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Guidance Law Development for Aeroassisted 
Transfer Vehicles Using Matched Asymptotic Expansions 

Anthony J. Calise* and Nahum Melamed** 

Georgia Institute of Technology, Atlanta, GA 30332 

Summary 

This report addresses and clarifies a number of issues related to the Matched 
Asymptotic Expansion (MAE) analysis of skip trajectories, or any class of problems that 
give rise to inner layers that are not associated directly with satisfying boundary conditions. 
The procedure for matching inner and outer solutions, and using the composite solution to 
satisfy boundary conditions is developed and rigorously followed to obtain a set of 
algebraic equations for the problem of inclination change with minimum energy loss. 

A detailed evaluation of the zeroth order guidance algorithm for aeroassisted orbit 
transfer is performed. It is shown that by exploiting the structure of the MAE solution 
procedure, the original problem, which requires the solution of a set of 20 implicit algebraic 
equations, can be reduced to a problem of 6 implicit equations in 6 unknowns. A solution 
that is near optimal, requires a minimum of computation, and thus can be implemented in 
real time and on-board the vehicle, has been obtained. Guidance law implementation 
entails treating the current state as a new initial state and repetitively solving the zeroth order 
MAE problem to obtain the feedback controls. 

Finally, a general procedure is developed for constructing a MAE solution up to 
first order, of the Hamilton-Jacobi-Bellman equation based on the method of 
characteristics. The development is valid for a class of perturbation problems whose 
solution exhibits two-time-scale behavior. A regular expansion for problems of this type is 
shown to be inappropriate since it is not valid over a narrow range of the independent 
variable. That is, it is not uniformly valid. Of particular interest here is the manner in 
which matching and boundary conditions are enforced when the expansion is carried out to 
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first order. Two cases are distinguished - one where the left boundary condition coincides 
with, or lies to the right of, the singular region, and another one where the left boundary 
condition lies to the left of the singular region. A simple example is used to illustrate the 
procedure where the obtained solution is uniformly valid to 0(e 2 ). The potential 
application of this procedure to aeroassisted plane change is also described and partially 
evaluated. 



Section I 


Introduction 


I.l Foreword 


Aeroassisted transfer problems concern the optimal maneuver of a space vehicle 
operating in vacuum around a planet with occasional passage through its atmosphere. As 
a rule, aerodynamic force is a dissipative force which has the effect of decreasing the 
vehicle's total energy. However, in terms of the fuel consumption, an aerodynamic 
maneuver can be inserted to achieve a transfer of the vehicle from an initial orbit to a 
destination orbit advantageously, compared to a purely propulsive exoatmospheric 
maneuver. Typical examples are a plane rotation of a space vehicle orbiting about a planet, 
and an aerobrake maneuver, in which the vehicle is to be transferred from a highly 
energetic hyperbolic orbit to a low energy elliptical or circular orbit, thereby ensuring its 
capture by the planet's gravitational field. Aeroassisted transfer problems are characterized 
by a maneuver phase in which the motion is dominated by atmospheric forces, and entry 
and exit phases where the motion is dominated by gravitational and inertial forces. In 
general, propulsive force can be used in any of the phases although in the so called 
aeroglide problem, it is not used in the aerodynamic part of the maneuver. The 
optimization of problems such as the plane change maneuver consists of guiding the vehicle 
using aerodynamic force, while minimizing the energy loss. Propulsive forces are used 
only in the exoatmospheric parts to deorbit and to compensate for the energy loss at the end 
of the plane change maneuver. 

An objective in any guidance study related to aeroassisted orbit transfer vehicles is 
the development of solutions that are implementable in real time and on-board the vehicle. 
Therefore, solutions that are both near optimal and that require a minimum of computation 
are of primary interest. The ideal solution from a computational point of view, is the one 
which reduces the solution of the governing differential equations to a set of algebraic 
equations, thus eliminating the need for multiple shooting or quadrature. This solution 
constitutes the basis for the optimal guidance algorithm. 
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1.2 Relationship to Earlier Results 

An extensive survey paper presented by Mease 1 gives the current status on the 
optimization of aeroassisted orbit transfer trajectories. In view of low cost transportation 
being a key to the utilization and exploration of space, important issues such as payload 
mass delivery capability and aerodynamic heating considerations are discussed. 
Aeroassisted transfer trajectories give rise to a difficult optimization problem from a 
guidance point of view. In general, numerical methods are required for an exact solution 
although approximation methods can be employed to obtain analytic solutions, or to reduce 
the solution to a quadrature. 


1.2.1 Numerical Optimization 

Examples of numerical solutions to optimal aeroassisted orbit transfer problems 
may be found in Ref.'s. 2-6, and in earlier works cited in these references. In Ref.’s. 2-5 
a family of problems were studied in the context of optimal aeroassisted orbit transfer, 
including the minimization of the time integral of the flight path angle squared. The 
solution to this latter problem results in nearly grazing trajectories that take place in an 
altitude range where viscous effects are expected. It is shown that the nearly grazing 
solution is a useful engineering compromise between energy requirements and aerodynamic 
heating requirements. Optimization subject to a hard constraint on heating rate was 
considered in Ref. 6, by reformulation as a parameter optimization problem. 


1.2.2 Analytical Studies 

Approximation methods can be employed to obtain analytic solutions, or to reduce 
the solution to a quadrature. However, it is difficult to precisely satisf y terminal co nstraints 
using these approximate for ms, because the nature of aeroassiste d transf ers is that the 
controls (typically lift and bank angle) are most effective while the vehicle is well within the 
atmosphere. Adjustments near the end of the maneuver to reduce terminal errors rapidly 
lead to control saturation/" 



Ref.'s. 7-1 1 typify the studies that have been performed on the problem of optimal 
aeroassisted orbit plane change, and that are directed towards obtaining analytical results. 
The authors of Ref.'s. 7-9 are able to integrate the state and costate equations by assuming 
that a quantity, known as Loh's term, M(h,V), is constant over the trajectory. M(h,V) 
represents the sum of gravitational and inertial forces, and is nearly zero for the entry phase 
of the maneuver, but unfortunately undergoes a rapid variation during the exit phase. In 
Ref. 10 a regular perturbation method is used, in which the perturbation parameter is 
identified as the ratio of the atmospheric scale height to the planet radius, and solutions are 
presented up to the first order in the perturbation parameter. The solution approach 
requires that quadratures be performed at each control update for the first order correction, 
essential to account for variations in M pear the end of the trajectory. The approach applied 
in Ref. 1 1 is similar to that given in Ref. 10 except for the use of an alternate independent 
variable. Unfortunately, large control effort can be observed near the end of the trajectory 
to satisfy terminal constraints. The interesting feature in Ref.'s. 10 and 1 1 is that the zeroth 
order solution corresponds to the constant M approximation in Ref.'s. 7-9, for which 
complete analytic results are available. 


1.2.3 MAE Studies 

In Ref.'s. 12 and 13, a MAE analysis is performed in which the perturbation 
parameter is the same as that used in Ref. 10. In Ref. 12 the expressions for the matching 
conditions are simplified over those obtained in earlier studies by using the inner solution 
alone to satisfy initial conditions. In Ref. 13 the state equations are integrated under the 
assumption of constant controls, and the expressions for the matching conditions are 
simplified by using the outer solution to satisfy initial conditions. The matching procedure 
taken in Ref. 12 is used to obtain an optimal lift control solution of an atmospheric entry 
problem, and in Ref. 13 to approximate an atmospheric skip trajectory for fixed controls. 
However, these approximate matching approaches are not recommended here for guidance 
law development since they are valid only when the initial condition lies either far outside 
or well inside the atmosphere. 

Singular perturbation methods for re-entry and aeroassisted transfer trajectories 
have, in more recent times, been explored in Ref.'s. 14-18. Ref.'s. 14-16 consider the re- 
entry problem. In Ref. 14 a singular perturbation parameter is artificially introduced on the 
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left hand sides of the altitude and flight path angle equations of motion. Several analytical 
guidance algorithms are derived for re-entry and evaluated by comparison to numerically 
obtained optimal solutions. While this approach appears useful for re-entry problems, it 
does not yield a satisfactory solution for aeroassisted transfer problems, because the 
boundary layer dynamics associated with satisfying the terminal constraints, is intractable. 
A suboptimal guidance algorithm is derived and evaluated in Ref. 15, which can be used in 
conjunction with the re-entry solutions in Ref. 14 to satisfy terminal constraints associated 
with aeroassisted transfer problems. In Ref. 16, a MAE analysis is performed in which the 
perturbation parameter is the same as that used in Ref. 10. The state equations are 
integrated under the assumption of constant controls, and the expressions for the matching 
conditions are simplified over those obtained in earlier studies by using the outer solution 
alone to satisfy initial conditions. Again, this approach is not recommended here for 
guidance law development, since it can only be used as an approximation when the initial 
condition lies far outside the atmosphere. 

Aeroassisted transfer problems are characterized by a maneuver phase, in which the 
motion is dominated by atmospheric forces, and entry and exit phases, where the motion is 
dominated by gravitational and inertial forces. The method of Matched Asymptotic 
Expansion (MAE) analysis is a mathematical realization of this intuitive description, and it 
is fundamentally different from re-entry problems in that the boundary conditions are given 
outside the region of singularity where the inner layer occurs. This type of problem gives 
rise to inner layers (regions where state and control variables can exhibit rapid variation) 
that correspond to intervals where the maneuver is dominated by aerodynamic forces. 
Thus the inner solution is associated with the aeroassisted portion of the maneuver. 

The application of MAE to aeroassisted transfer trajectories has been addressed in 
Ref.'s. 17 and 18. In Ref. 17, a general optimization problem is considered, and analytical 
results for the costate equations are presented. The problem is reduced to a set of constants 
of integration which are to be used to satisfy transversality conditions. Unfortunately, 
since the transversality conditions involve both states and costates, and since the state 
equations are not tractable in the inner layer, a multiple shooting method would have to be 
used to determine these parameters. Ref. 18 considers the development of an atmospheric 
guidance law for planar skip trajectories in which the controls are treated as constants to be 
updated at each guidance computation. This reference also identifies several 
inconsistencies that were encountered in the analysis, which are a consequence of an 
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incorrect approach in the application of MAE present both in Ref.'s. 17 and 18, and in 
earlier references taken from these papers. 


v 


1.2.4 MAE Analysis of the HJB Equation 

In Ref. 1 1 a regular perturbation analysis was performed for the aeroassisted 
inclination change problem and guided solutions were obtained in which the control 
solution was corrected to first order by regular expansion of the Hamilton-Jacobi-Bellman 
(HJB) equation. This work is similar to the results reported in Ref. 10, wherein the 
procedure was originally developed. The expansion entails quadrature along a 
characteristic curve, which is defined by the zeroth order Euler solution for the same 
problem. The quadrature must be repeated at each update of the control solution. A 
significant feature identified in this report is that the zeroth order solution from the regular 
expansion in Ref. 1 1 corresponds the zeroth order inner solution of the MAE formulation 
(but with modified boundary conditions), and that the zeroth order MAE composite 
solution represents a significant improvement over the zeroth order solution of Ref. 1 1 . 
This can be explained by the fact that the effect of gravitational and inertial forces are 
accounted for in the zeroth order MAE problem, whereas they are not in the zeroth order 
regular expansion solution where Loh's term is modeled as zero. The variations in this 
term are accounted for in the zeroth order outer problem of the MAE formulation. Of 
particular interest in this report is the manner in which matching and boundary conditions 
are enforced when the expansion of the HJB equation is carried out to first order. This is 
significantly more complex than for the case of regular expansion in Ref.'s. 10 and 11. In 
this report a general procedure for constructing a matched asymptotic expansion of the HJB 
equation, based on the method of characteristics, is developed for the first time. 

In Ref. 22, a uniformly valid power series solution to the HJB equation was 
obtained for a class of nonlinear, singularly perturbed systems, in which a small parameter 
appears on the left hand side of the equations of motion. The dynamics for this type of 
singular perturbation formulation is composed of slow and fast state variables, and the 
process of matching and forming a composition solution is considerably simplified by the 
fact that the left and right boundary conditions occur within the regions of singularity. For 
this class of singularly perturbed systems the slow and fast variables separately satisfy their 
respective boundary conditions in the zeroth order problem. In contrast, the MAE 
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formulation treated here is characterized by the fact that at any instant in the trajectory, the 
left boundary condition (which is the current state of the vehicle) may occur either to the 
left, within, or to the right of the region of singularity. For the problem of aeroassisted 
plane-change, this region reduces to a single point, corresponding to the lowest altitude 
point, in the limit as the perturbation parameter approaches zero. Moreover, it is not 
possible to identify separate slow and fast variables. For this class of systems, the problem 
is singularly perturbed due to the explicit dependence that the dynamics has on both t and 
t/e on the right hand side of the equations of motion, where t is the independent variable 
and e is the perturbation parameter. In the aeroassisted plane-change problem this 
parameter is closely approximated by the ratio of the atmospheric scale height to the radius 
of the Earth. 


1.3 Contributions 

This report addresses and clarifies a number of issues related to MAE analysis of 
atmospheric skip trajectories, or any class of problems that give rise to inner layers that are 
not associated directly with satisfying boundary conditions. A key point that has been 
overlooked in previous studies is that the inner solution is crucially involved in satisfying 
the boundary conditions as well as the outer solution. Other key issues that are either 
lacking, or improperly dealt with, in earlier studies are: (1) the need for left and right outer 
solutions, (2) the role that the inner solution plays in joining the discontinuities that occur 
between the outer solutions through the matching conditions, and (3) the need to properly 
select the reference altitude used in defining the independent variable of integration in the 
outer solution. The procedure for matching inner and outer solutions, and using the 
composite solution to satisfy boundary conditions, is rigorously followed in developing for 
the first time a complete algebraic solution to the problem of inclination change with 
minimum energy loss, uniformly valid to O(e). 

The motivation in the research reported here is that, by its nature, the aeroassisted 
transfer problem is better suited to analysis by singular perturbations rather than by regular 
perturbation analysis. In fact, it is shown that the regular perturbation problem in Ref.'s. 
10 and 1 1 is actually the inner part of a two time scale analysis based on singular 
perturbation theory. The interesting feature here is that Loh's term variations are now 
accounted for in the zeroth order outer solution (where the motion is Keplerian), for which 



an analytic solution can be obtained. The zeroth order inner solution corresponds to the 
solution presented in Ref. 9, but with important differences that pertain to the treatment of 
boundary conditions in a MAE approach. Moreover, it is shown that this approach is more 
accurate in satisfying terminal constraints than the regular perturbation approach in Ref.'s. 
10 and 11. 

The consequence of the above results is that the problem has been reduced to a set 
of algebraic equations, whose solution forms the basis for a feedback guidance algorithm. 
We have developed a zeroth order guidance algorithm, which is based on the MAE method, 
and is used to perform a detailed evaluation for the aeroassisted inclination change problem. 
Repeated solution of the algebraic equations resulting from the MAE analysis along the 
trajectory, treating each current state as an initial state, constitutes a feedback guidance 
algorithm. 

Finally, a general procedure is developed for constructing a MAE of the HJB 
equation based on the method of characteristics. Of particular interest here is the manner in 
which matching and boundary conditions are enforced when the expansion is carried out to 
first order, and how the characteristic curves are to be determined. It is shown that the 
MAE solution to the Euler system of equations associated with the zeroth order MAE 
problem provides the characteristic curves for the first order expansion of the HJB 
equation. Two cases are distinguished - one where the left boundary condition coincides 
with, or lies to the right of, the singular region and another one where the left boundary 
condition lies to the left of the singular region. A simple example is used to illustrate the 
procedure, where the obtained solution is uniformly valid to 0(e 2 ). The procedure's 
potential application to aeroassisted plane change is also described and partially evaluated. 

Ref.'s. 34-40 are all the publications that have resulted from this research. 


1.4 Outline of the Report 

Section 2 gives a general description of our approach in applying MAE to 
aeroassisted transfer trajectory analysis, and establishes some notation to be followed 
throughout the report. Section 3 uses the approach outlined in Section 2 in analyzing the 
orbit transfer problem of inclination change with minimum energy loss in the atmosphere. 
The analysis makes use of the problem formulation and analytic results in Ref. 14 for the 
outer solution. Then a transformation of variables is made to solve the inner problem, 
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using the analytical results presented in Ref. 9. This is followed by a description of the 
zeroth order MAE solution procedure, which makes use of a unique approach for solving 
the resulting algebraic equations. Numerical results are given to evaluate the resulting 
solution, which indicate a need for a first order correction. Section 4 gives the procedure 
for constructing a matched asymptotic expansion of the HJB equation to first order, based 
on the method of characteristics. It is also shown in Section 4 how to determine the 
characteristic curves of the expansion for two distinct cases, using a simple example to 
illustrate the procedure. This is followed by a detailed description of the procedure's 
potential application to aeroassisted plane change. Appendices A-C give details pertaining 
to intermediate derivations. 



Nomenclature 


a Integration Constant 

c Integration Constant 

C Aerodynamic Coefficient 

E Lift to Drag Ratio 

g Gravitational Acceleration 

h Normalized Altitude - Outer Independent Variable 
H Hamiltonian 

J Performance Index 

k Inner Integration Constants 

m Mass 

M Loh's Term 

P Return Function 

r Radius from Earth’s Center 

s Reference Area 

t Independent Variable, Time 

u Transformed Velocity - Outer Variable 

v, v Transformed Velocity - Inner Variable 

V Velocity 

w Transformed Altitude - Inner Variable 

x Dependent Variable 

Greek Letters 

P Inverse Scale Height 

8 Vertical Control Component 

e A Small Parameter 

4> Cross Range Angle 

0 Down Range Angle 

T| Stretched Normalized Altitude 
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Nomenclature (cont.) 

Y Heading 

y Flight Path Angle 

X Normalized Lift Coefficient 

|i Bank Angle 

p Atmospheric Density 

0 Horizontal Control Component 

X Independent Variable, Stretched Time 

Subscripts 

D Drag 

f Final Value 

1 Initial Value 

L Lift 

s Radius and Acceleration at Reference Altitude 
u Partial Derivative with Respect to u 

J Partial Derivative with Respect to y 

Y Partial Derivative with Respect to y 

Superscripts 

c Composite Variable 

i Inner Variable 

L Left Side 

o Outer Variable 

R Right Side 

* Corresponds to Maximum Lift to Drag Ratio 
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Section II 


Matched Asymptotic Expansion 
Analysis Procedure 


This Section clarifies a number of issues related to the MAE analysis of skip 
trajectories, or any class of problems that give rise to inner layers that are not associated 
directly with satisfying boundary conditions. In addition, a simple example is presented to 
^ illustrate the procedure. 


2.1 Conceptual Layout 

In any aeroassisted transfer problem, the vehicle passes through two distinct 
regions, in terms of the dominating forces (Fig's. 2.1 and 2.2). 




v 



Figure 2. 1 . Aeroassisted Orbit Transfer 
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Gravitational Forces Dominate 



Figure 2.2. Typical Skip Trajectory 

In the high altitude (outer) region, gravitational and inertial forces dominate, and the 
motion can be approximated as Keplerian, with the atmospheric effect considered as a 
perturbation. The low altitude (inner) region is dominated by aerodynamic forces, with the 
gravitational and the inertial effects considered as perturbations. It is crucial to note that the 
initial (approaching) and final (retreating) parts of the transfer trajectory are distinctly 
different in their trajectory parameters (when approximated as Keplerian arcs) as a 
consequence of what occurs in the aerodynamically dominated region. 

The method of MAE can be used as a mathematical realization of the above intuitive 
description 26 * 27 . It decomposes the total problem into two simpler subproblems (known as 
the inner and outer problems), which are appropriate for the separate regions of the total 
trajectory. The process of matching the solutions of these subproblems, and forming a 
composite solution, accounts for the perturbing effects as well, in a mathematically precise 
manner. For purposes of guidance law development, the objective is to obtain approximate 
analytical solutions in the outer and inner regions separately, and then to combine them to 
form a composite solution which is uniformly valid for the entire maneuver. 

To elaborate on this idea, consider the system of equations 

dx/dt = f(x , t, t/e, e) (2.1) 


The function f is assumed analytic with respect to its arguments in the region of interest, 
and in addition it is assumed to have the property that lim fix, t, t/e, e) exists for t * 0. 



The problem is singular in that the limit is not defined at t = 0. Conceptually, optimal 
control problem formulations may appear in this form after eliminating the control using the 
optimality condition, so that only state and costate variables appear in the equations of 
motion. This results in a two point boundary value problem, with left and right boundary 
conditions of the form tj) = 0 and ^Cxf, tf) = 0. Of particular interest here is the 
situation where tj < 0 < tf. 

The solution of Eq. (2.1) is sought in the form of an asymptotic series in e 
x°(t;e) = x°(t) + ex “(f) + e 2 x°(t) + (2.2) 


which is referred to as the "outer" expansion. The leading term in (2.2) is obtained as the 
solution of (2.1) for e=0. However, due to the singularity in f at t = 0, it is not a uniformly 
valid O(e) approximation of x(t, e). Note that for the situation tj < 0 < tf, two outer 
expansions are required, one for t < 0, and one for t > 0. 

Only zeroth order solutions are used in the analysis that follows in Chapters 2 and 
3, and to distinguish between the left and right zeroth order outer solutions, they are 
denoted by L Xg(t) and *r^(r), where the superscript o denotes outer, the subscript 0 

denotes zeroth order, and the superscripts L and R distinguish the left and right solutions. 
These solution segments are illustrated in Fig. 2.3, where it is important to note that in 
general there is a discontinuity at t=0. 
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To examine the solution in the neighborhood of t = 0, Eq. (2.1) is expressed in 
terms of a stretched independent variable r — tie 

dx/dx = ef(x, ex, x, e ) (2.3) 

The solution of Eq. (2.3) is sought in the form of an asymptotic series in e 

x'(r,e) = jf‘(r) + £xj(T) + e 2 x[(x) + ••• (2.4) 


which is referred to as the "inner" expansion. Again, only the leading term in Eq. (2.4) is 
considered in this Chapter and in Chapter 3. 

Matching the inner solution to the outer solution is performed separately for the left 
and right parts of the outer solution. To zeroth order in e, this is accomplished by 
enforcing the following relationships: 

xl(-oo) = L *°(0), 4(oo) = X°(0) (2.5) 

where x ' 0 (-°«), l Xq (0), x' 0 (°°) and *x% (0) are obtained by taking appropriate limits in their 
respective arguments x and t. The consequence of Eq. (2.5) on the limit behavior of x' 0 (x) 
is illustrated in Fig. 2.3, where the solution is shown superimposed with x£ ( t ) in terms of 
the original independent variable, t. The limit values in Eq. (2.5) also define the common 
parts of the inner and outer zeroth order solutions. 

At this stage a uniformly valid composite approximation can be constructed by the 
method of additive composition. The additive composition is obtained by taking the sum of 
the outer and inner solutions and subtracting the common part. In the left part of the 
trajectory, the composite solution is given by: 

X(';e) = %(0 + x' 0 (t/€) - L JCo(0), t< 0 (2.6) . 


and in the right part: 
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X(f;e) =X(t) + 4(t/e) -X( 0), t> 0 (2.7) 

From Eq's. (2.6) and (2.7) it is seen that at t=0 the composite solution takes the form: 

L x 0 c (0;e) =X(0;e) = *j(0) (2.8) 

Thus the composite solution is continuous at t=0. 

From the above discussion it is apparent that the inner solution (and the subsequent 
matching and forming of the composite solution) can be viewed as a process whereby the 
discontinuity between X(0 ar >d **o(0 at t= 0 is taken into account. In a skip trajectory, 

this discontinuity is caused by the change that occurs in the trajectory parameters during the 
osculating atmospheric portion of the maneuver. This point was not realized for example in 
Ref. 17 where a single outer solution was defined. The resulting composite solution is 
illustrated by the continuous bold line in Fig. 2.3. 

From Fig. 2.3 it is also apparent that the composite solution must be used to satisfy 
the boundary conditions. In Ref.’s. 16-18 only the outer solution was used to satisfy the 
boundary conditions at t=tj. In general, the contribution that x' 0 (t/e) makes to the 

boundary conditions depends on t; and tf. When the boundary conditions are far away 
from the region of influence of the inner solution, this contribution may be small. 
However, in an optimization problem it is well known that solutions exhibit large 
sensitivity to the boundary values of the costate variables. Also, in the context of 
developing a guidance algorithm, the initial condition should always be viewed as 
occurring anywhere along the trajectory. In Ref. 18 use of the outer solution alone to 
satisfy initial conditions led to discrepancies in the matching conditions which could not be 
completely resolved. 

In the analysis that follows, the convention of using 

t = h = (r - r,)/r, (2.9) 

is adopted, where r s is a reference radius. While h is not monotonic, it is apparent from 
Fig. 2.3 that this presents no conceptual difficulty in the outer solution, since separate 
constants of integration are defined for entry and exit. It is however important to transform 
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to a monotonic independent variable when integrating the inner dynamics. For this reason 
the transformation used in Ref. 9 is employed. The left and right composite solutions in 
Eq's. (2.6) and (2.7) are then distinguished by the sign of the flight path angle. 

Note that it is immediately obvious from Fig. 2.3 that t=h=0 corresponds to r s = 
fmin» where r m in is the minimum radius over the trajectory. Thus it is essential to use the 
condition 

rS(0;£) = yj( 0) = 0 (2.10) 


to determine r s , in contrast to the arbitrary selections made in Ref.’s. 17 and 18. 

The approach followed in the next chapter is first to express the inner and outer 
solutions in terms of constants of integration. The constants of integration for the left outer 
solution will be viewed as the variables used to ultimately satisfy the left and right 
boundary conditions. The constants of integration for the inner solution will be viewed as 
unknowns to be determined in terms of the left outer solution constants, using the left 
matching condition in Eq. (2.5). The constants of integration for the right outer solution 
are in turn evaluated in terms of the inner solution constants, using the right matching 
condition in Eq. (2.5). Finally, the following relationship, which holds when the j th 
component of x is constant in the outer solution, is to be noted 

x c Qj (t;e) = 4j(t/e) (2.11) 

and the following relationship which holds when the j th component of x is constant in the 
inner solution 


x c 0j (r,e) = x° 0j (t) (2.12) 

which follows directly from Eq's. (2.6) and (2.7). These relationships are used several 
times in the analysis that follows. 



2.2 A Simple Example 


2.2.1 Problem Formulation 

The exact dynamic is given by the following equation 

x = x + ue~" c / e, x(t t ) given, t f = 1, lul<l (2.13) 

where x is a scalar state variable and u is the control. Find a control u to maximize the 
performance index J=Xf. Using the maximum principle to evaluate the optimal control 
yields uOP^I and the exact solutions of the state and costate are: 

x( f,x(f.),f.) = [(1 + e)x(t i )e'" i + - «f ,/£ ]/( 1 + e) (2.14) 

A(r;jr(r,),r.) = e'~' (2.15) 


2.2.2 Outer Solution 

The zeroth order outer dynamic is obtained by taking the limit e 0 in Eq. (2.1 2) 

dx°Jdt = jc 0 ° (2.16) 

It is clear that this limit is not defined when t=0(e). This is precisely what makes this 
problem singularly perturbed and it should not be expected that the zeroth order outer 
solution alone will be a zeroth order uniformly valid approximation in t in the entire domain 
of interest. Specifically, it is not valid in the sub-region where t=0(e). We will elaborate 
on this more in Section 4. 

The zeroth order outer Hamiltonian is defined as 

Ho = Kx° 0 (2.17) 


and the zeroth order outer costate equation is given by 
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drjdt = -aH°/a* 0 ° = -a° 0 (2.1 8 ) 

Eq's. (2. 16) and (2. 1 8) can be integrated to obtain the zeroth order outer solution: 

x°o (0 = ae' (2.19) 

A°(0 = be- (2.20) 

where a and b are integration constants. 

2.2.3 Inner Solution 

In the inner problem a new stretched independent variable is defined as 
x = t/e (2.21) 

and the inner exact dynamics is obtained by using Eq. (2.21) in Eq. (2.13) 

dx' / dx = ex' + e~ T (2.22) 

The zeroth order inner dynamics is obtained by taking the limit e— >0 in Eq. (2.22) 

dx'Jdx = f- T (2.23) 

It is clear that Eq. (2.23) is not a zeroth order uniformly valid approximation in x to the 
exact dynamics in Eq. (2.22), because the exact dynamics is unstable, whereas the zeroth 
order dynamics is stable. Hence, for large x, the zeroth order inner solution is not a valid 
approximation to the exact solution. 

The zeroth order inner Hamiltonian is defined as 


h; = xy’ 


(2.24) 
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and the zeroth order inner costate equation is given by 

dX'Jdt = -9H’ /arj = 0 (2.25) 

Eq's. (2.23) and (2.25) can be integrated to give the zeroth order inner solution: 

jc'(r) = -e~ x + c (2.26) 

A’ 0 (t) = d (2.27) 

where c and d are integration constants. 


2.2.4 Zeroth Order Matching Condition 

The first step to remove the non uniformity from the zeroth order solution is to 
carry out the process of the matching the outer and inner solutions. To zeroth order the 
matching conditions entails equating the outer solution evaluated at small t (t=0) with the 
inner solution evaluated at large I ( r — » °°): 

a = x° 0 (0) = x' 0 (oo) = c (2.28) 

b = A° 0 (0) = A' 0 (~) = d (2.29) 


2.2.5 Zeroth Order Composite Solution 

The composite solution is constructed by the method of additive composition by 
taking the sum of the outer and inner solutions and subtracting the common part 

*0 (0 = *o°(0 + *o('/ £ ) - *o(0) (2.30) 


A c 0 (f) = K(t) + A ° 0 (t/e) - A° 0 (0) 


(2.31) 
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Making use of Eq's. (2.19),( 2.20) and (2.26- 2.29) in Eq's. (2.30) and (2.31) results in: 


*o(0 = ae' - e~' /c 

(2.32) 

Aq(0 = be~' 

(2.33) 


2.2.6 Boundary Conditions 

The composite solution is uniformly valid to zeroth order, hence it is used to 
enforce the boundary conditions. The initial condition on the state x is given in Eq. (2. 1 3). 
Using it in Eq. (2.32) evaluated at tj and solving for a yields 

a = e~' i [x(t i ) + e~ ,,u ] (2.34) 

Using this in Eq. (3.32) the zeroth order composite solution for the state variable becomes 
* 0 c (r; *(*.),*,) = elixir.) + e~' ,,e y - e~" c (2.35) 

The boundary condition on the costate is found from the definition of the performance 
index, defined as J=xf, as 

be~ l = Aq(1) = BJ/dx, = 1 (2.36) 

Solving for b and using it in Eq. (2.33), the zeroth order composite solution for X becomes 
A'tojr^),^) = e’-' (2.37) 

Note that since the zeroth order inner costate solution is constant, the composite solution is 
equal to the outer solution, and in this example it also equals the exact solution in Eq. 
(2.15). 

The individual zeroth order outer and inner solutions are: 
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V? 

x° 0 (t,x(tM) = MO + 

(2.38) 


A°o(?MO>0 = 

(2.39) 


xi(r>x(0,0 = -e~ T + [jc(0 + e~ T, ]e~ tl 

(2.40) 


A‘(r,x(0,0 = c 1 

(2.41) 


v Note that neither the inner or the outer zeroth order solutions for x, nor the inner zeroth 

order solution for X, are uniformly valid zeroth order approximations to the exact solution. 
For the costate this is obvious since the outer solution in Eq. (2.39), which is also the 
composite solution, is equal to the exact solution in Eq. (2.15). The inner solution for X in 
Eq. (2.41) is a zeroth order approximation to the exact solution only for the sub-region 
t=0(e). To show that the composite solution for x is a zeroth order uniformly valid 
approximation to the exact solution in Eq. (2.14), e is set to zero in Eq. (2.14) to get the 
zeroth order term in the expansion of the exact solution 

Um[x(r;jr(r, ),/,)] = <?“'■[ Jt(r,) + e~' ,lc ]e' - e~" c (2.42) 

which is precisely the zeroth order composite solution for x given in Eq. (2.35). 

V 

2.2.7 Location of the Initial Condition 

v The singularity in this example is located at t=0. Fig 2.4 gives the solution for the 

state x in the case £=0. 1 , tj=0 and x(tj)=0 that is, the initial condition is given at the 
singularity. Eq. (2.30) indicates that at t=0 the composite solution equals the inner 
solution. Hence the inner solution and the composite solution, but not the outer solution, 
v satisfy the initial condition as is apparent in Fig. 2.4. It is also apparent that only the 

composite solution is a uniformly valid approximation to the exact solution and that the 
matching condition is satisfied. 

In realistic situations the initial conditions, or the boundary conditions, do not in 
^ general occur at the singularity (see the discussion after Eq. (2.8)). Fig 2.5 gives the 
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solution for the state x in the case e=0. 1 , tj=0.2 and x(tj)=0.5, that is the initial condition 
and the singularity do not coincide. It is evident that neither the inner solution nor the outer 
solution satisfy the initial condition. The composite solution does satisfy the initial 
condition and it is the only solution that uniformly approximates the exact solution, as 
before. Also note that the region of singularity is best approximated by the inner solution 
alone, as before. This situation basically represents either the left or the right side in Fig. 
2.3 for aeroassisted maneuvers. It clarifies why two sets of matching condition are needed 
in aeroassisted maneuvers, and the role that the inner solution plays in joining the 
discontinuity of the outer solution to give a continuos composite solution at the singularity. 



Figure 2.4 Exact, Composite, Outer and Inner State Solutions 
for initial Condition x(0.0)=0.0 



Figure 2.5 Exact, Composite, Outer and Inner State 

Solutions for initial Condition x(0.2)=0.5 





Section III 


Inclination Change With Minimum Energy Loss 


In this Section the procedure for matching inner and outer solutions, and using the 
zero order composite solution to satisfy boundary conditions is rigorously followed in 
developing a complete algebraic solution to the problem of inclination change with 
minimum energy loss. Repeated solution of these algebraic equations along the trajectory, 
treating each current state as an initial state, constitutes a zero order feedback guidance 
algorithm. 


3.1 Problem Formulation 

The three dimensional point mass equations of motion for a lifting vehicle over a 
spherical non-rotating planet are given by: 


dr / dt - Vsiny (3.1) 

dd / dt = V cosy cos y/rcosQ (3.2) 

d<f>/dt = Vcosysin^/r (3.3) 

dV / dt = - psC D V 2 / 2m - gsiny (3.4) 

Vdy/dt = psC L V 2 cos// / 2m - (g - V 2 /r)cosy (3.5) 

Vdy/dt - psCjV^ sin p /2m cosy - V 2 cosy cos y tan <p /r (3.6) 
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where the coordinates system is defined in Fig. 3. 1 



A Newtonian gravitational field has the form 

g(r) = g/jr 2 (3.7) 


where the subscript s denotes a reference radius, defined to be the minimum trajectory 
radius. An exponential atmosphere model is used: 

P(r) = p t e~* lr ~ r '\ p = MH t (3.8) 


where H s is the scale height The lift and drag coefficients are assumed to be of the form 
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C L — C L A 

(3.9) 


C D = C*(l + A 2 )/ 2 

(3.10) 


where the constants C* and C*> are the lift and drag coefficients corresponding to the 
maximum lift to drag ratio, and X is the normalized lift coefficient. 

Defining the following dimensionless quantities: 


h = (r - r,)/r. 

(3.11) 

« = V 2 /8S> 

(3.12) 

B - C *P* S / 2m P 

(3.13) 

* * , * 


E = C L /C„ 

(3.14) 

e = UPr , = HJr, 

(3.15) 


and using h instead of t as the independent variable, the state equations can be written as 
follows: 



dd/dh = cos^coty/(l + /i)cos0 

(3.16) 

V* 

d(f>/dh = sin y^cot //(l + h) 

(3.17) 


du/ dh = - Bu( 1 + A 2 )e~ hfc / eE* sin y - 2 / (1 + hf 

(3.18) 


dy/dh = B\cosne~ h,£ /esiny + [1/(1 + h) - \/u(l + h) 2 ]coty (3.19) 
dy/ / dh = BXsinfj£~ h ^ £ It sin y cosy - cos^tan0coty/(l + h) (3.20) 
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The objective is to achieve an inclination change with minimum fuel consumption. 
In Ref. 7 it is shown that for short duration maneuvers, the change in the cross range angle 
$ is small, and the change in inclination is closely approximated by the change in the 

heading. Moreover, fuel consumption is nearly minimized by minimizing the energy loss 
in the atmospheric phase of the maneuver. Furthermore, in Ref. 19 it is shown that 
although the actual inclination change depends on the initial inclination, the starting point of 
the maneuver can be timed so as to obtain the maximum inclination change that is achieved 
when the initial inclination is zero. A consequence of this result is that the initial plane can 
be taken as the plane of reference, which will be refered to as the equatorial plane. Under 
these assumptions, 0(0) and <|>(0) can be set to zero without loss of generality, and the 
heading angle \|/ is used to approximate the change in the inclination angle. Thus© and <j) 
become ignorable coordinates, and the equations of motion reduce to: 

du/dh = - Bu{ 1 + A 2 )e~ hlc / eE* sin y - 2/(1 + hf (3.21) 

dy/dh = BA cos ne~ hlc / e sin y + [1/(1 + h) - 1/m(1 + hf]coiy (3.22) 

d\j/ / dh = BX sin fie~ hlc /csinycosy (3.23) 

The parameter e is the ratio of the atmospheric scale height to the minimum trajectory 

radius. In general, r s is nearly equal to the planet's radius, and, for the purpose of 
calculating e to form the composite solution, it will be treated as such. For Earth, the value 
of e is of order 10‘3. The controls are the normalized lift coefficient A and the bank angle 

Using the definition of u given in Eq. (3.12), the objective of minimizing energy 
loss can be equivalendy expressed as: 

max(J}, J = u f (3.24) 

The Hamiltonian function associated with the state Eq's. (3.21-3.23) has the form: 


H = [-Bu(l + X 2 )e~ h,c / eE* sin y - 2/(1 + h) 2 ]P u 



+ {BA cos iie~ hlc / esin y + [1/(1 + h) - l/u(l + h) 2 ]coty)P r 
+ [BA sin \ie~* u / esin ycos y]P w (3.25) 

where p u » Py and py are the associated costate variables. Assuming the controls X and |i 
are not beyond their limits, their optimal values are obtained as a function of the state and 
costate variables by solving the optimality conditions Hx.=0 and H^=0. The resulting 
expressions are: 


A = E*(P y cos ^ + P ¥ sin n / cos y) / 2 uP u 

(3.26) 

tan fj. = P y / P y cos y 

(3.27) 


3.2 Zero Order Outer Solution 17 

The zero order equations for the outer problem can be obtained by simply taking the 


limit as £ approaches zero on the right hand side of Eq's. (3.21-3.23): 

du°Jdh = -2/(1 + hf (3.28) 

dfjdh = [1/(1 + h) - l/nj(l + hf~\cotyl (3.29) 

dyTJdh = 0 (3.30) 

The general solution for the outer system to zero order in e is given by: 

= 2[c, + 1/(1 + A)] (3.31) 

cosy 0 °(A) = c 2 /(l+ A)Vuo°(A) (3.32) 


Vo(h) = c 3 


(3.33) 
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where ci , C 2 and C 3 are constants of integration. The adjoint equations are given by: 

dPJdh = -H x (3.34) 

where x is any of the state variables. In the outer region these equations to zero order in e 
are: 

dP° u Jdh = -/^cotyJ/Kd + h)f (3.35) 

dP;jdh = P° r0 [ 1/(1 + h) - l/u 0 °(l + A) 2 ]/sin 2 y° (3.36) 

dP; Q /dh = 0 (3.37) 

The solution to this system using Eq's. (3.31-3.33) is: 

P° u0 (h) = -a 2 /2uo° + a, (3.38) 

P; 0 (A) = a 2 tany 0 ° (3.39) 

1 

P°r 0 (h) = a 3 (3.40) 

where aj , a 2 , and 83 are constants of integration. Eq's. (3.31-3.33) provide the exact 
solution in the outer region and are the integrals of Keplerian motion that express 
conservation of energy and conservation of angular momentum. 


33 Zero Order Inner Solution 9 

In the inner region, where the aerodynamic force is dominant, a new stretched 
altitude is defined as: 



29 


*7 * h < E (3.41) 

However, in order to integrate the inner layer equations the transformations used in Ref. 9 
are also adopted. These have the additional feature of transforming the independent 
variable from Tj to a monotonic independent variable, y. Analytic solution also requires a 
small flight path angle approximation. 

cosy=l, sin y = y (3.42) 

The following transformations are defined 20 : 


w = Be~" 
v = E* ln(l / g t r t u) 

The controls are also transformed to vertical and horizontal components: 

8 = Acos/i 
o = A sin fi 

Invoking the above transformations, Eq's. (3.21-3.23) become: 

dy/dy = [8 + eM(e)]/o 
dw/dyr = -y/o 


(3.43) 

(3.44) 

(3.45) 

(3.46) 

(3.47) 

(3.48) 


dv/dy/ = [1 + 8 2 + tr 2 + eG(e)]/o 


(3.49) 
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where em(e) is Loh's term which accounts for the gravitational and inertial effects on the 
motion. 


Af(e) = {1 - e v,E 'g t r t /[ 1 + eln(B/w)]}/H{l + eln(B/w)] (3.50) 

and 


G(e) = e VIE g,r,£*sin2y/w[l + fin (B/w)f (3.51) 

The approximations used in Ref. 9 to integrate the state and costate equations were 
that eM(e) is constant over the trajectory, and that the term eG(e) is negligible. We note 
here that setting e = 0 in Eq's. (3.47) and (3.49), to obtain the zero order inner solution, 
corresponds to the approximations in Ref. 9 with em(e) = 0. Actually, EM is nearly zero 
during entry but undergoes a sharp variation during the exit phase. This is the main 
shortcoming to the approximation in Ref. 9. The interesting feature in this analysis is that 
the zero order inner solution corresponds to the M = 0 solution in Ref. 9, but that M 
variations are accounted for in the outer solution. It should also be pointed out that while 
the integrated solution bears a close resemblance to that in Ref. 9, it is totally different in 
the method of evaluating the constants of integration. 

The transformed zero order state equations in the inner region will be obtained by 
letting e=0 in Eq’s. (3.47) and (3.49). The zero order state equations for y' 0 and v’ 

become: 


dy\ / dyi = 5/ o 

(3.52) 

dvl/dx// ' 0 = [1 + 5 2 + o 1 ]/ a 

(3.53) 


Using Eq. (3.44), the performance index Eq. (3.24), expressed in the transformed 
variables, takes the form: 


-V. /£* 

— a * 


J = £ 




max{J}, 


(3.54) 



Note that uf in Eq. (3.24) and vf in Eq. (3.54) represent the composite solution of the 
transformed final velocity to zero order, and not the inner or outer solutions alone. The 
zero order Hamiltonian function in the inner region is given by: 


Hj, = P'o r 5/o - Pl w y / a + />*„( 1 + 5 2 + CT 2 )/ CT (3.55) 

and is constant since it does not depend on the independent variable y/' 0 . The control 
functions are obtained from the optimality conditions H' 05 = 0 and = 0. They are 
given by: 

5 = -P' 0r /2Pi (3.56) 

o 2 = i-Plf -P'ZyKMPl + 1 (3.57) 


The corresponding transformed zero order costate equations in the inner region are obtained 
using Eq. (3.34): 


dPo r /dy' Q = P\Jc 

(3.58) 

dPIJdri = 0 

(3.59) 

dPi/dy ’ = 0 

(3.60) 


In Appendix A it is shown that HJj = 2oP' v0 , and a consequence of this result and Eq. 
(3.60) is that o is constant too. 

At this stage the state Eq's. (3.48,3.52,3.53) and the costate Eq's. (3.58-3.60) can 
be integrated with respect to y/' 0 to result in: 
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m'o(V'o) = ttiVo /6 - ^ / 2 - *3V4 + * 4 ]/<r (3.62) 

Vo(V'o) = (<T+ 1 /<r)yr* + oKyj*, -^) 3 ]/3^, + k 5 (3.63) 

and: 

Plyivl) = (Pi I a) Vo + c (3.64) 

f*L = constant (3.65) 

Pi = constant (3.66) 


where k 3 , k 4 , k.s and c, P' Qw , P' 0v are the constants of integration of the state and costate 
equations respectively. The constants kj and k 2 are defined in terms of these constants as: 

k, = Pi Ha 1 Pi, k, « -c/laPi (3.67) 


3.4 Matching Conditions 

The method of additive composition (as described in Section 2), is used to combine 
the zero order inner and outer solutions into a single, uniformly valid approximation. The 
additive composition is obtained by taking the sum of the solutions in the different regions 
and subtracting the common part. Matching implies agreement between the outer solution 
for small values of h (h— >0), and the inner solution for large values of r) ( r/ — > <»). Since 

the inner solution lies between two discontinuous outer solutions, matching is done 
separately in the left and right parts of the transfer trajectory, and two sets of matching 
conditions result. Each set of conditions involves matching of both the state and costate 
variables. 



3.4.1 State Matching Conditions 


The solution for the states in the outer region is given by Eq's. (3.31-3.33). 


Taking the limit h— >0 yields: 


X(0) = 2(c, L + 1) 

(3.68) 

cos[ L y°(0)] = c\ / [2(Cj L + 1)] 1/2 

(3.69) 

Vo(0) = c\ 

(3.70) 


where superscript L denotes the constants of integration for the left outer solution. 

The solutions for the transformed inner variables are expressed with y/' 0 as the 

independent variable. Conceptually, it is possible to perform an inverse transformation to 
express the inner solution in the original variables with t] as the independent variable. 
However, the inverse transformation can be bypassed by first recognizing from Eq's. 
(3.43, 3.44) that 


L wi(oo) = 0 (3.71) 

X(~) = E*ln[l/g /s X(oo)] = E* ln[l / 2g t r t (c, L + 1)] (3.72) 

where the second relationship in Eq. (3.72) follows from Eq. (3.68) and enforcement of 
the matching condition X (0) = X(°°)- Applying the matching condition to y' 0 and ys' 0 

gives: 


7o(~) = l 7o°(0) = cos -1 {cj / [2(c, L + 1)] ,/2 } 

(3.73) 

¥o(°°) =X(0) = c 3 l 

(3.74) 


Using Eq's. (3.71-3.74) to evaluate Eq's. (3.61-3.63) at rj = °o yields: 
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cos" , {c 2 l /[2(c 1 l + 1)] 1/2 J = -k^cf / 2 + ^c 3 l + k 3 (3.75) 

0 = kjC^/6 - kff 1 2 - k 3 c 3 + k A (3.76) 

£’ ln[l / 2g,r,(c, L + 1)] = (a + 1 1 o)c\ + oKc, 1 *,-*,) 3 ]/^ + * 5 (3.77) 


Eq's. (3.75-3.77) can be used to evaluate k 3 , k 4 and k s in terms of the constants of 
integration for the state variables in the left side outer solution. Recall that k, and k 2 have 
been defined in Eq. (3.67). 

An exactly symmetric set of equations results from the matching conditions for the 
state variables on the right side: 

cos"" 1 {cj / [2(c* + 1)J ,/2 } = - k,cf / 2 + itjC* + k 3 (3.78) 

0 = V* 3 / 6 ~ h c ? I 1 - k 3 c] + k 4 (3.79) 

£*ln[l/2g/ t (cf + 1)] = (a+ \/o)c\ + o[(clk l -k 2 ) i ]/3k i + k s (3.80) 

Eq's. (3.78-3.80) relate c*, c“ and c 3 to the constants of integration for the state variables 
in the inner solution. 


3.4.2 Costate Matching Conditions 

The left and right matching conditions for the costate variables are defined below. 
The inner solution for the costate variables, corresponding to the transformed state 
variables, is given in Eq's. (3.64-3.66) in terms of the constants of integration c, P' 0w and 
P' 0v . To perform the matching with the outer costate solutions given in Eq's. (3.38-3.40), 
it is necessary to first find the corresponding expressions for P' 0u and P' 0¥ for the inner 

solution. In effect, this amounts to transforming the inner solution back to the original 
problem variables. 



The value of any costate variable at time t can be interpreted as the sensitivity of J to 
perturbations in the corresponding state variable at tune t (Ref. 17). Thus it can be written: 


P u = di/du ; P v = dJ/dv (3.81) 

From Eq. (3.44) it follows that: 

dv/du = -E* /u (3.82) 

Combining Eq's. (3.81) and (3.82) it becomes apparent that 

PL = PivOv/34 = -E*PL/K (3.83) 

To determine P' 0v ,\t is noted that \ff' 0 is the independent variable in the transformed 

inner problem. Therefore, making use of the Hamilton-Jacobi-Bellman equation 17 results 
in: 

Pi,* = dJ/av'i = -h; (3.84) 

In Appendix A it is shown that H' 0 = 2 crP' 0v , and thus: 

P' 0¥ = ~2oP' 0v (3.85) 

Using Eq's. (3.83) and (3.85), the costate matching process is ready to be carried 
out . Taking the limit h-+0 in Eq's. (3.38-3.40) and using Eq's. (3.68-3.70) the left side 
matching conditions results: 

L PL( 0) = -a\IMc\ + \)+ a, L = L P' 0 J~) (3.86) 

L P°o r ( 0) = a l tan{cos '{c2 / [2(c, L + 1)] ,/2 }} = L /»’ r (oo) (3.87) 
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^oV 0 ) = = l /V~) ( 3.88) 

Taking the limit V-*°° in Eq's. (3.64, 3.83, 3.85) and making use of Eq's. (3.68, 3.70, 
3.74) and of the state matching condition L u° (0) = X («*>) leads to: 

L /i(~) = -£>' v /2( Cj L + 1) (3.89) 

l ^o r (°°) = (f*L / ^3 + c (3.90) 

^orH = -2cPL (3.91) 

Substituting Eq's. (3.89-3.91) in Eq's. (3.86-3.88) yields: 

-a^!A(c\ + ]) + a\ = -E'Pl!2(c\ + 1) (3.92) 

tan {cos -1 {^2 /[2(c, L + 1)] J/2 )} = ( PL/o)c\ + c (3.93) 

a\ = -2 cPl (3.94) 

Finally, in Appendix A it is shown that a can be written as: 

a = 1/(1 + X + 2*,* 3 ) ,/2 (3.95) 

Eq's. (3.92-3.94) relate c, and P' 0v to the constants of integration for the left outer 
costate solution. 

A precisely symmetrical set of equations results from the matching conditions for 
the costate variables on the right side: 

-a} / 4( c* + 1) + af = -E’Pil2(c\ + \) (3.96) 



tan{cos -1 {c* / [2(c" + 1)] ,/2 }) = (/>' w / a)c\ + c 


(3.97) 


a\ = -2cK 


Ov 


(3.98) 


Eq s. (3.96-3.98) relate af , a\, and a* to the constants of integration for the inner costate 
solution. 


3.5 Zeroth Order Composite Solution 

The composite solution is constructed by the method of additive composition as 
outlined in Section 2. The form of the composite solution for the state and the costate 
variables is: 


*§(/!,£) = xqW + xq ( h/e) - jcq(O) (3.99) 

where now e = H s /f s , f s is the planet radius (see the comment following Eq. (3.23)) and x 
is any of the state or costate variables. The solution for the outer state variables is given in 
Eq’s. (3.31-3.33) and for the inner state variables in Eq’s. (3.61-3.63). The outer left 
solution, evaluated at h=0, is given in Eq's. (3.68-3.70). Using Eq. (3.44), the left zeroth 
order state composite solution is given by: 

X(M) = 2/(1 + h) + *f v ° ( * /t)/£ ‘ / g s r t - 2 (3.100) 

Yo(h,e) = cos"'{c 2 L /(l + /t)[2(c, L + 1/(1 + h))] m ) 

+ {-^[¥o(h/ e)f /2 + k^ylih/e) + * 3 } 

- cos -1 {Cj / [2(Cj L + 1)] 1/2 } (3.101) 

Vo(M) = Vo(h/£) (3.102) 

Note that the composite solution for y is simply the inner solution since the outer solution 
is constant (see the general comment at the end of Section 2). Since no outer constants of 
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integration appear explicitly in the left composite solutions of u and of y, the left and right 
composite solutions are identical. The right zeroth order composite solution for y is given 
by: 


Yo(h,E) = cos H {c* / (1 + A)[2(c* + 1/(1 + h))) ]n ) 

+ {-^[ylih/e)] 2 / 2 + /fcjV'o (*/*) + * 3 } 

- cos -1 {Cj / [2(c* + 1)] 1/2 } (3.103) 

By Eq. (3.43), w^, as a function of h/e is given as: 

w i(h/e) = Be' h,c (3.104) 

and Eq. (3.62) is used to solve for yr' 0 (h / e) in terms of the altitude h and the constants of 
integration. Eq. (3.63) is then used to evaluate v' 0 (h/ e). Thus Eq's. (3.100-3.103) 

together with Eq's. (3.62, 3.63, 3.104), provide the composite solution for the state 
variables as a f u nction of h. 

The form of the zeroth order composite solution of the costate variables is given in 
Eq. (3.99) and it is constructed the same way the composite solution for the states was 
obtained. The solution for the outer costate variables is given in Eq's. (3.38-3.40) and for 
the inner costate variables in Eq's. (3.64-3.66). The outer left solution evaluated at h=0 is 
given in Eq’s. (3.86-3.88). Using Eq's. (3.44) and (3.83, 3.85), the left composite 
solution is given by: 

L P c 0u (h,£) = ~a\! 4[Cj L + 1/(1 + h)] 

-PIE' l[e- v ' Wc)IE ' / g,r,] + oj /4(c[ + 1) (3.105) 

l P 0 c r (/t,£) = tan{cos" ] {c 2 /(I + h)2 m [c\ + 1/(1 + A)] 1/2 }} 

+ (Pl/°)vUh/e)+ c 

- tan{cos _1 {cj / 2 1/2 [c, L + 1]‘ /2 }} (3.106) 


X ¥ (h,e) = a\ 


(3.107) 



The values of v' Q (h/ e) and \f/' 0 (h/ e) are found from Eq's. (3.62, 3.63) and (3.104). 
Eq's. (3.40) and (3.85) for P 0¥ show that it is constant in the outer and inner regions. 
Hence the composite solution for P 0¥ is also constant and is simply the outer costate 
constant of integration. On the right side, the composite solution takes the form: 

'PL(h,E) = + 1/(1 + h)] 

-P'lX l[e-^ (hlc)IE ' /g.r,] + aj /4(c; + 1) (3.108) 

'P c 07 (h,e) = aj tan{cos _1 {Cj /(I + h)2 m [ C] “ + 1/(1 + /i)] 1/2 }} 

+ {KJ<y)v[(,hle)+ c 

- tan{cos -1 {cj / 2 1/2 [c* + l] l/2 }} (3.109) 

'P c 0¥ (h,e) ■ oj (3.110) 


From Eq's. (3.94) and (3.98) it is seen that a\ = a\, thus the zeroth order composite 
solution for P 0¥ is constant all along the transfer trajectory. Eq's. (3.105-3.1 10) provide 

the zeroth order composite solution for the costate variables as a function of h. 

The reference radius, r s , is also treated as an unknown parameter in the problem. 
As discussed in Section 2, it is the distance to the lowest point of the transfer trajectory at 
which h=0. Using the fact that at this point the composite solution for y is zero, provides 
the relationship needed to evaluate r s . The composite solution for y is given in Eq's. 
(3. 101) and (3. 103). At h=0 it becomes the inner solution evaluated at h=0, 

Yo(0,e) = 7o(0,e) = Xo(0) = o (3.111) 

The solution for y' 0 is given in Eq. (3.61). Equating it to 0 gives a relationship for the 
value of y/' 0 corresponding to h=0, 
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and the solution for y/' G (0) is: 


V'o(O) = *,/*,[ 1 ± (1 + 2^ /^) ,/2 ] 


(3.113) 


This value of yo(0) is used in Eq. (3.62) to evaluate w' 0 at h=0 


wj(0) = [*,y’( 0) 3 /6 - ^^o( 0) 2 /2 - *^(0) + * 4 ]/cr (3.114) 


Finally, use of Eq's. (3.1 14) and (3.13) in Eq. (3.104) for h=0 yields: 

Wp(0) = Cj) s s/2mP (3.115) 

from which p s (and thus r s ) can be found in terms of the constants of integration of the 
problem. 


3.6 Enforcing Boundary Conditions 

To complete the solution of the problem, the initial conditions and transversality 
conditions at the end point must be satisfied. As discussed in Section 2, this is performed 
by using the composite solution, and can be thought of as the process by which the 
constants of integration for the left outer solution are evaluated. 

Assume that the initial conditions u;, Yi and yj are given, along with the initial 
radius r,. For the purpose of exposition, assume that Yf, yf and the final radius r f are given. 
The corresponding transversality condition is: 

P e Ah f ) = 1 (3.116) 


where 
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V 


V 






^ = in ~ r,)/r tl h f = ( r f - r,)/r, (3.117) 

Using u,, Y„ \|fj, Yf> Yf and Eq. (3.1 16), the following equations result from enforcing the 


boundary conditions: 

u, = 2/(1 + h,) + <r vi( * l/£)/E ’ /g t r t - 2 (3.118) 

Yi = cos -1 [cl / (1 + tome] + 1/(1 + h,))] 112 ) 

+ hWty / e)} 2 / 2 + kM(Me) + k 3 ) 

- cos -1 {Cj / [2(C| L + 1)) ,/2 } (3.119) 

¥i = Vo(h,/e) (3.120) 

V, = Yo(*//£) (3-121) 

1 = -a\ / 4[Cj R + 1/(1 + h f )] 

-PIE' /[ e - v ° ih ' lc)IE ‘ /g t rj + a\ /4« + 1) (3.122) 

r f = cos -1 {Cj / (1 + Ay)[2(c* + 1/(1 + h f ))] in ) 

+ / e)] 2 / 2 + k 2 y' 0 (h f / e) + k 3 ) 

- cos -'Icj/Pfc; + 1)] I/2 1 (3.123) 


where hj and hf are defined in Eq. (3.1 17). To summarize, Eq's. (3.75-3.80, 3.92-3.98, 
3.115, 3.118-3.123) constitute a set of 20 equations for the 18 unknown constants of 
integration c\> c\, c 3 \ k 3> * 4 , * 5 , c,\ c*, c 3 , a, L , a£, tf 3 , c, Pq v , a*, a', a], and 

for the parameters c and r s . 




42 


3.7 Zeroth Order MAE Numerical Solution 

3.7.1 Solution Procedure 

In this section, by exploiting the structure of the MAE solution procedure, it is 
shown how to simplify the original problem by further reducing it to a set of 6 implicit 
equations in 6 unknowns. The unknowns are the common parts of the inner and outer 
solution (which are equal to one another). 

The iteration procedure involves repeated solution of the inner and outer problems 
using the common parts as artificial boundary conditions. The common parts are adjusted 
in the iteration process until the actual boundary conditions are satisfied by the composite 
solution. The matching conditions are enforced at each iteration by simply equating the 
inner and outer common parts. We also exploit the fact that heading is constant in the outer 
solution, which permits enforcement of the actual boundary conditions on heading using 
the inner solution alone. This allows enforcement of two of the six boundary conditions at 
each iteration, and further reduces the problem to four equations in four unknowns, for 
which a Newton method is used to obtain a solution. Representative zeroth order MAE 
approximate solutions were calculated to exhibit the resulting inner, outer and composite 
solutions. 

The vehicle used in this study is the Maneuverable Research Re-entry Vehicle 
(MRRV). Its aerodynamic and mass characteristics are presented in Ref. 24. The 
maximum lift to drag ratio is 2.362, the lift coefficient corresponding to (L/D) m ax is 
Cl=. 1512, the reference area is s=l 1.69 m 2 and the mass is m=4898.7 kg. The constants 
needed for the exponential atmospheric density function are obtained by fitting to the 
standard atmosphere densities at 30 and 60 km. The resulting value of the scale height is 
H s = 1/^=7625.4 m. 

A Newton method was first tried to solve for the 20 unknowns. This approach was 
not successful due to the complexity of the relationships. An alternative approach was then 
derived wherein the number of coupled equations was reduced to 6 by defining the 
unknowns to be the common parts of the separate inner and outer solutions. The equations 
that determine these unknowns are the original boundary conditions enforced on the 
composite solutions. The basic idea is to first use the common parts as artificial boundary 
conditions to evaluate the constants of integration in Eq's. (3.31-3.33, 3.38-3.40) and 
(3.61-3.66), starting with an initial guess. Then a Newton method is used to iterate on the 



common parts until the original boundary conditions are satisfied by the composite 
solution. 

This procedure was slightly modified to capitalize on the structure of the MAE 
solution for this problem. Since y/o 1S constant in the outer solution, it follows from Eq. 
(3.102) and the matching conditions that y/ c 0 (the zeroth order composite solution for y) is 
simply the zeroth order inner solution ( y/' Q ). This allows the boundary conditions on y to 
be enforced using only the inner solution at each stage of the iteration process. Nested 
Newton iteration algorithms where implemented to determine the 6 unknowns using the 
boundary conditions for the composite solution. The inner iteration procedure corresponds 
to solving the inner problem with the boundary conditions on y enforced. The outer 
iteration process is used to enforce the four remaining boundary conditions on the 
composite solution. 

The first step in the procedure is to calculate the 6 inner solution integration 
constants k 3 , k 4 , k.s, c, P' 0w , and P' 0v in Eq's (3.61-3.66) using a guess for the common 
part values of on the left and for the common part values of y' 0 , P' 0v on the right 

(the constants kj and k 2 are given in terms of the other constants in Eq. (3.67)). Note that 
an initial guess for the common parts of w' 0 is not needed since the common part 

corresponds to rj — > Hence it follows from Eq. (3.43) that the common part values for 
Wp are zero on both the left and right sides. This particular arrangement concerning which 

of the common part values are treated as unknowns was chosen to agree with the actual 
boundary conditions for the original problem, which greatly simplifies the problem of 
forming an initial guess. The procedure of forming an initial guess and the equations for 
evaluating the constants of integration are given in Appendix B. 

Next, an inner loop Newton search is performed on the left and right common part 
values of y/{, (using only the inner solution) so that the boundary conditions = Vi 

and y„(wy) = V/ are satisfied, where wj and Wf correspond to r = rj and r = rf in Eq's. 

(3.43) and (3.1 17). This is done by solving the cubic equation in Eq. (3.62) while taking 
into account that w' 0 (yf) for the range of y of interest must be positive for the solution to 

have physical meaning. The inner solution procedure also determines r s and a (see 
Appendix B). 

Fig. 4. 1 presents an example of a converged solution of this equation for yi=0°, 
yf=20°. The resulting left and right common parts of yr\ are -1.10° and 20.75° 
respectively, which correspond to w' 0 = 0 ( rj -» «>) in this figure. The values of at y = 
0° and y = 20° map to z, = zf = 60 km altitude above sea level via Eq. (3.43). These 
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altitudes where specified as part of the boundary conditions. Thus the boundary conditions 
Yi( r i) = 0° and Y f (r f ) = 20° have been met. 



-5 0 5 10 15 20 25 

V (deg) 


Figure 3.2. Example of Converged Inner Newton Iteration 
Process for w as a Function of y 


The left and right outer solution integration constants in Eq's. (3.31-3.33) and 
(3.38-3.40) are obtained by enforcing the matching conditions in Eq's. (3.75-3.80) and 
(3.92-3.94, 3.96-3.98). This amounts to equating the common parts for the outer solution 
(Eq’s. (3.31-3.33, 3.38-3.40) evaluated at h=0) with the common parts from the inner 
solution (Eq's. (3.61-3.66) evaluated at the left and right common part values of y' 0 
corresponding to *4=0 in Fig. 3.2). Since some of the outer solution variables are 
different from the inner solution variables, it is necessary to employ the transformations 
in Eq's. (3.43, 3.44) and (3.83). Note that C 3 and a 3 are not needed since the outer 
solutions for y/'o and Po r arc constant. Hence the zeroth order composite solution for 

these variables is the inner solution alone. The zeroth order composite solution is then 
evaluated to determine if the boundary conditions are satisfied: 
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V 


S' 


- «, = 0, Xo (A ) - Yi = 0, 

PLify) -1 = 0, rS^) - Yf = 0 (3.124) 

The error in these equations is used in an outer loop Newton iteration to iterate on the 
remaining four unknown common parts. These are the left common part values for u and 
y, and the right common part values for P„ and y. A summary of this procedure is given in 
Fig 3.3 
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Initial Guess of Common Part Values 




Evaluation of Inner Integration Constants ^ 

1 

Newton Iterations for the left and 
right Common Parts of V 

i 

Evaluation of Outer Integration Constants 

I 

Enforcing Boundary Conditions Using Composite Solution 



Update Estimate of 
Common Part Values 


Exit When Accuracy Requirement is Met 


v 


w 


Figure 3.3. Summary of the Numerical Solution Procedure 
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3.7.2 Open Loop Numerical Results 

Fig's. 3.4-3.7 present the zeroth order MAE converged solution for the boundary 
conditions: Zj = Zf = 60 km (initial and final altitudes above sea level), Vj=7851 m/sec, y,=~ 
1.35°, Yf=l 0° and A\jr = 20°. The corresponding control time histories are given in Fig's. 
3.8 and 3.9 and the altitude history in Fig. 3.10. The value of the expansion parameter 
used here is e = l/(3r s = 1 1.87e-4. In Fig's. 3.4-3. 7 the inner solution is shown for the 
complete range of iy between its left and right common part values, which correspond to 
7 /— >oo ( w' 0 =0 in Fig. 3.2). This was done to illustrate the fact that the matching 

conditions in Eq's. (3.75-3.80) and (3.92-3.94, 3.96-3.98) are satisfied on the left and 
right portions of the solution (see also Fig. 2.3) 

These results clearly indicate that the composite solutions for y and py are 
significantly different from the inner solution. In particular, the major variation in py is due 
to the outer solution, which in effect amounts to a correction for the large variation in Loh's 
term during the exit phase. We note that the outer solution plays an important role in 
forming a uniformly valid zeroth order approximation to the exact solution. The 
normalized lift control X is always near 1 .0, corresponding to flight at near maximum lift to 
drag ratio. The bank angle |i is always near 90° indicating that most of the aerodynamic 
force is utilized in performing the turn. 




y(deg) 

Figure 3.5. Zero Order MAE Inner, Outer and 
Composite Velocity Solution 
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V (deg) 

Figure 3. 10. Zero Order MAE Altitude History 


3.8 Zeroth Order Guided Solution 

3.8.1 Solution Procedure 

Zeroth order closed loop guided solutions are obtained by using the optimal control 
expressions given in Eq's (3.26) and (3.27). These expressions involve both the states 
and the costates, thus knowledge of both states and costates is needed to evaluate the 
controls. Assuming that the states are available for feedback, only estimates of the costates 
are required at each control computation along the trajectory. Feedback implementation 
entails treating the current state (from the simulation) at each control update as a new initial 
state, and calculating the costate values corresponding to the same time instant. The 
estimate for these costates (to zeroth order) are obtained by repetitively solving the zero 
order MAE problem. 

In the first step, an initial guess and boundary conditions are supplied to initiate the 
procedure of obtaining a zero order MAE converged solution (Appendix B). Next, the 
costate expressions in Eq's. (3.38-3.40) and (3.64-3.66) can be evaluated as a function of 




the corresponding independent variables and used in Eq's. (3.105-3.107) to construct the 
composite costate expressions. These are in turn used in Eq’s (3.26) and (3.27) to 
compute the controls. When a predetermined time increment has been reached, the current 
states are used as initial conditions for the next MAE calculation. It follows that the initial 
guess is available in every step of the zero order MAE calculation after the first 

For control computation between MAE solution updates, the integration constants 
from the last update are used. The transformations defined in Eq’s. (3.1 1, 3.12, 3.43) and 
(3.44) are used to transform the simulated dimensional variables into the inner and outer 
dependent and independent variables. The transformed altitude h is used in Eq’s. (3.31, 
3.32, 3.38) and (3.39) to compute the left and right outer costates. The heading y is used 
in Eq. (3.64) to evaluate P' 0v , and Eq's. (3.105-107) are used to calculate the composite 
costates. The composite costates and the current y and u (from Eq. (3.12)) are used in 

Eq's. (3.26) and (3.27) to evaluate the optimal controls between the time instants where the 
MAE solution is updated. 

During the exit phase, the left outer solution is discarded, and matching is required 
only between the right outer solution and the inner solution. In this case, the constants of 
integration for the inner solution are viewed as free parameters used to satisfy the boundary 
conditions. 


3.8.2 Numerical Results 

Fig's 3.11 through 3.14 present a comparison between the optimal solution 
(obtained using a multiple shooting method 2 -* 5 ) and the zero order guided solution for 
A\|/=20°. The corresponding control time histories are given in Fig's. 3.15 and 3.16. 
Loh's term (corresponding to the optimal solution) is given in Fig. 3.17 and the altitude in 
Fig. 3.18. The time increment between guided solution updates is 5 seconds, with the 
control updated at every integration step following the procedure described above. These 
results indicate that the guided solutions and the optimal solutions are in a very good 
agreement throughout the trajectory. The error that does exist in some of the guided 
solution variables, for example the y and the altitude (z) solutions, indicate the need for a 
first order correction. A general procedure for MAE expansion of the HJB equation to first 
order is developed in Section 4, along with a description of its potential application to aero- 
assisted orbit transfer. 
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Fig's 3. 13 and 3. 15 clearly indicate how the variation in Loh's term near the end of 
the trajectory is partially accounted for in the zero order guided solution. Specifically, the 
normalized lift coefficient X, given in Fig. 3.15 does not saturate in the exit phase, but 
reduces to near zero. This is due to the fact that the correction in Py from the outer solution 
is too large (see also Fig. 3.6). Fig. 3. 1 9 compares the velocity histones near the end of 
the trajectory. The optimal and guided solution values for terminal velocities are 6751 
nVsec and 6736 m/sec respectively. The difference of the guided value from the optimal one 
in 15 m/sec. 




























Section IV 


Matched Asymptotic Expansion 
of the Hamilton-Jacobi-Bellman Equation 


In this Section a general procedure for constructing a matched asymptotic expansion 
of the Hamilton-Jacobi-Bellman equation, based on the method of characteristics, is 
developed. The development is valid for a class of perturbation problems whose solution 
exhibits two-time-scale behavior. A regular expansion for problems of this type is shown 
to be inappropriate since it is not valid over a narrow range of the independent variable. 
That is, it is not uniformly valid. Of particular interest here is the manner in which 
matching and boundary conditions are enforced when the expansion is carried out to first 
order. Two cases are distinguished - one where the left boundary condition coincides with 
or lies to the right of the singular region and another one where the left boundary condition 
lies to the left of the singular region. 


4.1 Singularly Perturbed Hamilton-Jacobi-Bellman Equation 

Ref.'s. 10 and 1 1 have formulated the aeroassisted plane change maneuver as a 
regular perturbation problem and employed a regular expansion of the optimal return 
function (P) to first order to obtain an approximate guidance law for that problem. 
However, as shown in Section 3, the optimal aeroassisted plane change maneuver is 
actually a singularly perturbed problem. For the solution of such a problem to be 
uniformly valid, inner and outer expansions of P and a matching procedure have to be 
performed so that a uniformly valid composite solution for the optimal return function can 
be constructed. This Section treats the expansion to first order for a class of singularly 
perturbed problems that are characterized by the presence of t and t/e in the right hand sides 
of the differential equations. 
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We assume that the dynamic system in Eq. (2.1) has the following form 
dx/dt = f(x,u,t) + g(x,u,t/ £)/ £, x(f,) given (4.1) 

where x e /?" is the state vector, u e R m is the control vector, e is a small parameter and t is 
the independent variable. The functions f and g are assumed analytic with respect to their 
arguments in the region of interest and, in addition, g is assumed to have the property that 

lim g(x,u,t / e) / e‘ =0, t* 0, / = !,... (4.2) 


The singular characteristic of this type of problem arises at t=0, where the above property is 
not satisfied. To investigate the behavior in the region of singularity near t=0, a stretched 
independent variable t is defined as x=t Jt and Eq. (4. 1) in terms of x is given by 

dx/dr = ef(x,\i,ez) + g(x,u,r), x(r.) given (4.3) 


where the function f is assumed to have the property that 

lim ef(x, u,£T) = 0 
£-►0 


(4.4) 


Note that both systems in Eq's. (4.1) and (4.3) represent the exact dynamics. Eq. (4.1) is 
called the outer system and all the variables associated with it will be denoted by the 
superscript o. Eq. (4.3) is called the inner system, and all the variables associated with it 
will be denoted by the superscript i. 

The assumption in Eq. (4.2) regarding the form of g was chosen because the 
aeroassisted plane change problem satisfies this property. It is typically satisfied when g in 
Eq. (4.1) has the form g = e~" c h(x,u). Also, the division by £ in Eq. (4.1) insures that 
the zeroth order inner dynamics in Eq. (4.3) are not zero. As will become evident, this 
specialization is not essential to the development in this section, and the main results are not 
restricted to it Generalizations (and specializations) to other assumed forms for both f and 
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g are straightforward extensions of the main ideas to be presented below. The only 
requirement is that the zeroth order subproblems remain well defined. 

Note that when the transformation to Eq. (4.3) is performed, it appears to result in a 
regular perturbation problem, when in fact the original dynamics in Eq. (4.1) are singular. 
This is precisely what was done in Ref.’s. 10 and 11. It will be shown later, in the 
application section, that transformation of our equations that begin in the form of Eq (4. 1 ) 
to the variables used in Ref. 1 1 results in a set of equations in the form of Eq. (4.3). This 
means that the inner expansion alone was used in these earlier studies to satisfy the 
boundary conditions, and that the expansions are valid only within the region of 
singularity. In an MAE expansion, the solution is dominated by the inner expansion within 
this region, because the outer expansion nearly cancels with the common part of the 
solution. Therefore, the differences will be most apparent outside this region, and for 
problems involving skipping trajectories through the atmosphere this will be most evident 
where Loh s term undergoes as large variation. That is, whenever the orbital forces 
dominate. These forces are ignored in the zeroth order regular expansion solution, and are 
accounted for in the zeroth order MAE solution. This point was illustrated in the numerical 
results in Section 3. 

The optimization problem is to find u(x,t) that minimizes J = <p{Xf \ subject to 
the dynamic constraints in Eq. (4.1) and the terminal constraint \y(x f ) = 0, where 
x f = x(t f ) and tf is the final time. The outer HJB equation is 

P ° = ~^ H = -W" + ^'/e) (4.5) 


where U is a class of continuous bounded controls, / opI s f(x,u or *,t), 
8 ° P ‘ s 8(x,t / c.m^'Cm)), and w opl (.*:,/ > I ,f) is given by the optimality condition H u = 0, 
assuming that H m is positive definite. P a {x 0 (t l ),t i ) is the optimal return function defined 
as the optimal value of the performance index for an optimal path starting at x(t l ) and r, 
and satisfying the terminal constraint. In terms of the stretched independent variable x the 
inner HJB equation is given by 


/>; = -min H = -ftaT + g°*) 


(4.6) 
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where / >, (r , (r,),r i ) is the same optimal return function as defined for Eq. (4.5), but 
expressed in terms of the inner variables. It is important to note that in both formulations 
the exact solutions for the return functions are identical and depend on the same initial 
condition x°(t i ) = x' ( r, ) s x^), where Xj=tj/e. This is not true for the corresponding inner 

and outer expansions of the solution, since the boundary conditions is satisfied with the 
composite solution. Therefore, only the composite solution should be compared with the 
exact solution. However, the inner and outer expansion solution terms may be viewed as 
being dependent on x(tj) and tj through the matching conditions and the boundary 
conditions enforced on the composite solution. This dependence will be made explicit in 
the following development 

Consider a power series expansion in e of the return functions for both the outer or 
inner formulations: 


P°(x(0,t„£) = X F J WM*' 

;» o 
/'- 0 


(4.7) 

(4.8) 


By substituting these expansions into the outer and inner expressions of the optimal control 
(derived from the optimality condition) and expanding in a power series in e, the outer and 
inner expansions of the optimal control in a feedback form are obtained as: 


u oop , ( x ( 0,/>;,0 = X ^ U )* ; = u ° 0 ( t Wt , U „£) 

y»o 


(4.9) 


u iap '(j c(r, ),/»;, o = £«;(*(*, xoe' = u iop ‘(x(t,),t„£) 


(4.10) 


y»o 


The details regarding this expansion for a regular perturbation problem are given in Ref. 
11, and they apply directly for the separate expansions in Eq's. (4.9) and (4.10). The 
zeroth order terms u£ and u' 0 are the optimal controls for the zeroth order outer and inner 

problems that are obtained by setting e=0 in Eq's. (4.1) and (4.3). If analytic solutions are 
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available for the zeroth order problem, then the higher order solutions are determined by the 
expansions of the outer and inner HJB equations: 


P: = I P°,e j = -(lPl'e i )(lf°(xJ)e i ) 

j »0 j »0 


(4.11) 


P\ = I Pj T e j = lg)(x,r)e i ) 

1-0 j -0 >« 0 j —0 


(4.12) 


where the expansions for f and g come from substituting the series in Eq's (4.9) and (4. 10) 
for u in Eq's. (4.1) and (4.3). The reader is again refered to Ref. 1 1 for the details, while 
noting the exception that the dependence of f in Eq. (4.3) on e enters both through u and 
ex. Equating like powers of e in Eq's. (4.1 1) and (4.12) leads to outer and inner first order 
linear partial differential equations for P° and P) 


p° + p° f° 
r j t ^ r jxJ 0 

= R](x(tM’PU’ 

no J 

•••»* o 

j=l,... 

(4.13) 

P 1 + P 1 p 1 

r jx ^ r jx6o 


P) 

•I 1 0 / » 

j=l,... 

(4.14) 


where the forcing terms R° and /?’ are functions of the lower terms in the solution. This 

procedure was carried out in Ref. 11 for a regular expansion, and is identical in form for 
the inner and outer expansions defined here with obvious accounting for differences in 
notation, and the exception noted above regarding the dependence of f on ex in the inner 
expansion. In particular, the j=l expressions are: 

K = 0 , R ; = -PiJoix, 0) (4.15) 


after taking into account the optimality condition for the zeroth order problem. Note that 
R° is zero as a consequence of the assumption in Eq. (4.2), which also accounts for the 
fact that an expansion of g is not required in Eq. (4.1 1). In the sequel it is not assumed that 
R° = 0 so as to allow for more general assumptions regarding the form for g. 
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Singular perturbation expansions differ from regular expansions in that solutions to 
any order of the outer and inner problems are not uniformly valid approximations of the 
exact solution. Uniformly valid approximations are constructed from the outer and inner 
solutions by the MAE method. This process involves matching of the inner and outer 
solutions and enforcing the boundary conditions on an additive composite solution made up 
of the outer solution plus the inner solution minus the common part. See Sections 2 and 3 
for details. 

The partial differential equations in Eq's. (4.13) and (4.14) may be solved by the 
method of characteristics (Ref. 18). The characteristic curves for any order term of 
P° and P' are defined by the zeroth order differential equations: 

k = roc«oV) ( 4 . 16 ). 

dx'Jdx = g'(x^tX) (4.17) 

whose solutions are denoted by x 0 o (r;x(r,),f,) and jc;(r,x(r,),r,) respectively. It is shown 
in Appendix C that the boundary conditions for Eq's. (4.16) and (4.17) are also defined by 
the MAE method. Therefore the characteristic curves for the inner and outer expansions are 
the inner and outer solution components of a zeroth order MAE solution. Note that to any 
order, the boundary conditions are satisfied by the composite solution and not individually 
by the inner and outer solution components. The same is true for the characteristic curves 
which satisfy Eq's. (4.16) and (4.17). However, both the inner and outer characteristic 
curves depend on x(tj) and tj through the matching conditions and the boundary conditions. 
See Sections 2 and 3 for details. Note that the expansions P° and P] in Eq's. (4.13) and 

(4.14) are useful only at the initial time t 0 = er 0 , because the characteristic curves are 
computed for a specified x(tj) and tj. 


4.1.1 First Order Solution for the Case tj£0 

In order to express the matching condition, it is necessary to express the outer 
solution of Eq. (4. 13) for the interval 0 £ t £ tf, and the inner solution of Eq. (4. 14) for the 
interval 0 < x <, We also assume in this section that the initial time satisfies r, £ 0, and 


V 

v— . - 
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will subsequently generalize the result for arbitrary tj. The solutions for P° 
expressed as: 

and P\ are 

- 

i 

= \R:(r,x(tM)dt + p;(0,x(tM) 
0 

(4.18) 


T 

Pl(r,x(tM) = \RU^ x (tM)dr + PliO-MtM) 
0 

(4.19) 

V 



: 

A necessary condition in MAE analysis is that the zeroth order solution at t=0 serves as a 
stable equilibrium point for the zeroth order inner solution. This means that the 
characteristic curve for the inner expansion must approach a well defined limit as 
Hence it is assumed that/^'Cr;*^),^) reaches a well defined limit as t— >=» which will be 
denoted by /?, (*(/, ),/,). Thus, in order to express P\ for large values of x, it is necessary 
to rewrite Eq. (4.19) in the following form: 

: 'fsF 

T 

Pl(r,x(tM) = j[R;(r,x(tM)-RndT+ tR? + P;(0; *(/,), 0 
0 

(4.20) 

: -w 

We are now ready to perform the matching procedure between the outer and the inner 
solutions of the return function by enforcing the following rule: 

V' 

/ >0 (smallt) = P' (large 

(4.21) 


where the dependency on *(/, ),/, is omitted to save notation. Retaining terms to first order 
in t and e yields 

: 

W) + Rl.m + <(0) = /»’( -) + £[/>;(0) + /[*;(?) -/?;** jdr + 

0 

- 


(4.22) 
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where Eq. (4.20) is used to express P' for large x. Equating like powers in t and e, the 


following matching conditions are obtained: 


e°: P 0 °( 0) = P' 0 (oo) 

(4.23) 

e': PU 0) = P\( 0) + JlRl(t) - R\~]dz 

0 

(4.24) 

Rf = P° 0I ( 0) = -F 0 V°(0) 

(4.25) 


Next the boundary conditions are applied to the composite solution expressed to 
first order. The composite solution is constructed by adding the outer and inner solutions 
(expressed to first order) and subtracting their common part. The common part is the left 
hand side of Eq (4.22). Evaluation of the resulting first order composite solution at the 
final time yields the following boundary condition 

W = W + e[P?(t f ) + - roor, + £P,°(0)]= Kx,)(4.26) 


where the zeroth order composite return function is given by 


W = ^o° iff) + P' 0 (T f ) - p° 0 ( 0) 


(4.27) 


From Eq's. (4.18) and (4.20) the first order terms evaluated at tf and tf are related to their 
respective values at t = x = 0 by: 


P°(t f ) = jR°(t)dt + P,°( 0) 
0 


(4.28) 


p l (*» 


|[/?;(x) - Ri~]dT+ Rfz f + />;(0) 


(4.29) 



Using Eq's. (4.25) and (4.27-4.29) in Eq. (4.26), and recalling that z f = t { / e yields: 


c°: Pld f ) = <P(x f ) (4.30) 

e 1 : P/(0) = -jR°(t)dt - /[/?;( r) - Rf]dz (4.31) 

0 0 

Eq's. (4.31) and (4.24) are used to evaluate the constants /^(O) and P\( 0) in Eq's. (4.18) 
and (4.20). Finally, using these results in Eq. (4.26) evaluated for an arbitrary initial time 
tj such that 0 < ti < tf results in 


= P'oixVM) + £{-lR:(r,x(tM)dt 

T / 

-/[^(r.^O.r.) - Rl~]dz) 


(4.32) 


Partial differentiation of this expression with respect to x(q) provides an approximation for 
the costate variables to first order, which, when used in the optimality condition, results in 
an approximation for the feedback control to first order. 11 


P^WM.e) = P^Jx(t,),t,) + £{-J R^Jt,x(tM)dt 

+ R°(t / \x(t i ),t i )dt / / 3jc(/,) -}[^ ) (r,x(r i ),0 - R;~Jdz 

Tf 

+ [RUz f ;x(tM) ~ RH^f/Mt,)) (4.33) 


4.1.2 A Simple Example 

Here, the simple example that was presented in Section 2.2 is followed up. The 
problem formulation is 
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x = x + ue~" c /e, Jc(f.) given, t f = 1, I«I<1 (4.34) 

where x is a scalar state variable and u is the control. Find a control u to maximize the 
performance index J=Xf. To repeat the main results, the exact solutions of the state and 


costate are: 

x( nxitiXti) = [(1 + e)x(t i )e'~' i + - e~" c ]/(l + e) (4.35) 

A (t f x(tM) = e‘~‘ (4.36) 

and the zeroth order outer and inner solutions are: 

x° 0 (f,x(t,),t,) = W<i) + e-‘ i,c ]e'-- (4.37) 

A° 0 (r;x(O,O = e'~' (4.38) 

4 (r;x(r ( ),r,) = -e~ x + [x(r,) + e' Xi ]e~'‘ (4.39) 

A' 0 (r,jc (f.),f,) = e' (4.40) 


In this case the zeroth order composite solution for X equals the outer solution since the 
inner solution is constant. Comparison of Eq's. (4.36) and (4.38) shows that it also 
equals to the exact solution. Hence, no correction terms for X to first and higher orders are 
expected. 

Using Eq. (4.15) it is found that 

R?(t) = 0 (4.41) 

Ri(r) = -PU t)x' q (t) = -e l [-e~ x + e~'‘(x(t,) + e~ ,JC )] (4.42) 

Evaluating R\ ( r) in Eq. (4.42) at t — » yields 
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w 


*r = + e~ ,Jc y\ (4.43) 

Note from Eq's. (4.37,4.38) and (4.43) that the matching condition in Eq. (4.25) is 
automatically satisfied. This will always be the case when the outer and inner zeroth order 
MAE solutions are used to define the characteristic curves (see Appendix Q. Using Eq's. 
(4.41-4.43) in Eq. (4.32) gives 


v? 


Me 

Pi(x(tM>£) = Po(t,)~£j [Rl(t)- RI~Vt 

tile 

Me 

= P 0 e (O -ee' je~ T dt = P 0 c (r,) + ee l [e~ Uc - e~' i,c ] (4.44) 

Ule 

Note that the first order term in Eq. (4.44) does not depend on x(t t ), hence there is no 
correction to first order for the costate function as expected. The exact return function is 
the final value of x from Eq. (4.35), with t replaced by tf =1. 

P(x(tM>£) = [(1 + £)Jc(li>‘‘* + - e~ llc ]/(\ + e) (4.45) 

To show that Eq. (4.44) is a uniformly valid approximation to 0(e), the exact 
solution in Eq. (4.45) is expanded to first order 

JWiH.*) = PoixitM) * £[e~ Uc - (4.46) 

The difference between the first order terms in Eq’s. (4.44) and (4.46) is 

J(e,r,) = ee l (l-e~‘ i )e~' i,c + ee~ Uc (l-e l ) = eE{e y t t ) (4.47) 

The range of interest for tj is 0 < tj ^1 . The region of singularity corresponds to q=0(e), or 
to the range 0 £ tj < ke, where k is some constant. Outside this region ti=0(l). We now 
investigate the size of £(£,/,) both inside and outside the region of singularity: 


v- 
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for tj=0(e): 

limE(£,r)/ £ = ke'~ k <<*> 
£-►0 

(4.48) 

for tj=0(l): 

lim£(£,f.)/£ = 0 
£~>0 

(4.49) 


Thus, £(e,/,)=0(e) when tj is O(e) close to zero and o(e) when tj is 0(1). Hence, A(e,ti) 
in Eq. (4.47) is 0(e 2 ) uniformly in tj. It follows that the solution in Eq. (4.44) is a 
uniformly valid first order approximation to the exact solution in Eq. (4.46). When the 
regular expansion of Ref. 1 1 is applied to this problem with x = t/e as the independent 
variable (this transforms the dynamics to the form of Eq. (4.3)), then the first order 
solution is not valid in the region outside the neighborhood of t=0. 


4.1.3 First Order Solution for the Case tj<0 

In this case the region of singularity near t=0 occurs between the initial time ti<0 
and the final time tf>0. As detailed in Sections 2 and 3, two outer expansions are required 
for this situation, one for the interval t < 0 and the other corresponding to t > 0. The outer 
expansions are in general discontinuous at t = 0, but the composite solution is always 
continuous 11 . The inner expansion must be considered for -«> < x < oo in order that 
matching may be performed on both sides. This is illustrated in Fig. 4.1. 



Figure 4. 1 . Inner and Outer Expansions of P in the Left and Right Regions 
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The equivalent expressions to Eq’s. (4. 1 8) and (4.20) on the left side for t<0 are: 



(4.50) 

(4.51) 


Note that, in general, the inner forcing function has different limit values on the left 
and right sides as denoted by the L and R superscripts. Matching on the left side is done in 
a fashion similar to that on the right side and the equivalent expression to Eq. (4.21) is 

P ° (small negative t) = P' (large negative r)| T (4.52) 



Retaining terms to first order in t and e gives 


V 

7> 0 °(0) + L />°,(0)r + £ l />;( 0) = 

p 'o(-°°) + etf(0) - }[/?;( T)- L /?r]rfT + 7?;“r] 

— M 

(4.53) 


Equating like powers in t and e yields: 


V' 

o 

£ 

0 

ll 

o'* 

✓*— V 

1 

8 

'w' 

(4.54) 


0 

e 1 : 7?(0) = p;c 0) - J[/?;(T) - 

— * m 

(4.55) 


t': L /?r =^ 0 ,( 0 ) = - L P°oMT(0) 

(4.56) 


Since the inner characteristic curve is continuous at t=0, so is P\( 0), and therefore it 
satisfies the right side expression in Eq. (4.31). Using Eq. (4.31) in Eq's. (4.51) and 
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(4.55), and then Eq. (4.55) in Eq. (4.50) gives the inner and outer solutions of the return 
function on the left side. The composite solution on the left side is constructed the same 
way Eq. (4.26) was constructed and, after some manipulation, it takes the form 

Plih) = + e{-jR°(t)dt + />i(0) - J[/?j ( t) - L Rrm (4.57) 

Finally, using Eq. (4.31) in Eq. (4.57) yields 

'/ 

= P c 0 (x(tM ) + ehjR^Ma^dt 

0 

- J[/f 1 i (r,x(f,),r i ) - L /?ridr - J[/?;(r;x(r < .),r i ) - */?r]dr} (4.58) 

T t 0 

where satisfies the right side matching condition in Eq. (4.25). Differentiating Eq. 
(4.58) with respect to x(t t ) yields 


'/ 

P L< X ('M) = PI (x(0,t,)+ 

h 

0 

- R;{t f \x{tM)dt f /dx(t,) - f[/?; (r ;x(t,\t,) -K ¥t 

J *<U *0.) 

r, 

x l 

- f[/?; (v,x(tM) -'Ri~ ]dz 

' *w 

- IMh)) (4.59) 


which, when used in the optimality condition, results in an approximation for the feedback 
control to first order. 

As noted earlier, the first order corrections in Eq’s. (4.33) and (4.59) are only 
useful at to, and therefore the quadratures must be repeated at each control update. 



4.2 Application to Aeroassisted Plane-Change 

In Section 3 analytic zeroth order solutions were obtained for the aeroassisted 
plane-change problem. In this problem, the singularity occurs at the lowest altitude point, 
hence the result in Eq. (4.59) must be used for initial conditions when they correspond to 
the entry phase of the maneuver. As shown in Section 3, the reduced three dimensional 
point mass equations of motion (3.21-3.23) for a lifting vehicle over a spherical non- 
rotating planet have been written in the form of Eq. (4.1). Note that these equations satisfy 
the condition in Eq. (4.2) and that h plays the role of the independent variable, which is 
identified as t in Eq. (4.1). It is a monotonic decreasing variable to the left of the 
singularity (entry phase of the trajectory) and a monotonic increasing variable to the right of 
the singularity (exit phase of the trajectory). Therefore t in Eq. (4.3) is h/e, which is 
denoted by T]=h/e in Eq. (3.41). 

The control expressions in Eq’s. (3.26) and (3.27) depend on P y , P Y and P„. 
These can be obtained to first order using Eq. (4.59). From the first of Eq. (4. 1 5) R® = 0, 

hence only the inner quadrature is needed. From the second of Eq. (4.15) and Eq’s. 
(3.21)-(3.23) 

=^ou(-2) + Poy[0“l/“i)/tanyj] (4.60) 

^i°° in e Q- (4.59) can be obtained by using the matching conditions 
^ou(°°) = ^om( 0)> ^oy(°°) = ^oy(O), Wo(°°) = “o (0) and the expressions in Eq's. (3.31), 
(3.38) and (3.39). The result is simply 

R *r(~)=2VV l /C(oo ) = 2 VS (4.61) 

Since tany] appears in the denominator of Eq. (4.60) and yj(0) = 0 [see Eq. 
(3.111)], a singularity occurs in the integrand of Eq. (4.59) for this problem at T]=0. 
However, the singularity is removable by simply transforming the variable of integration 
0i) to the independent variable for the inner problem (y). Using Eq's. (3.23), (3.41) and 
(3.43) it follows that 


dr\ = sin y cos ydy / ow 


(4.62) 
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which removes the singularity in Eq. (4.59) when the variable of integration is changed 
from T| to \\f. 

Let yr s denote Vq ( rj = 0). This quantity is known from the zeroth order solution. 
Then, the specialization of Eq. (4.59) to the AOTV problem becomes: 



W f i i 

- / [-%**(*) - d\ff +%*i(¥0 - *Rhi /dVi) (4.63) 

y s <^0 'Vf 


Vs J 

?!, (>,) = rL «i) +«(- 1 H-M*{(v0 - L K'n\,iv 

r ‘ r ‘ o 

Vf J . j , | 

- I HWlW - R /T]| dVf/dYi) (4.64) 

^ ovv 0 ow 0 1 V/ 

rf u w = n u W +*{-T[At^(r) - L *nMv' 

0 ^0 

- rfr -tA^i(r) -"«r]]| (4.65) 

Vs <™ 0 O»o l Vf 1 

where a small flight path angle approximation y * sin ycos y is employed for simplicity of 
notation. 

The zeroth order composite terms in Eq's. (4.63-4.65) are known from the zeroth 
order solution as developed in Section 3. The expression for Rj(y) follows from Eq. 
(4.60), Eq's. (3.44) and (3.83) that relate to and Pl u to P\ v , and Eq's. (3.61), 
(3.63), (3.64) and (3.67) that relate yj(vO> Vq(^), P' 0r (Y) and P' ov to y and the 

constants of integration. R'(y/) is evaluated along the inner characteristic curve, which 
corresponds to the zeroth order inner solution in Section 3. Hence, differentiation or Rj 
with respect to x(q) implies that dxl / dx{t{) must be computed. The dependency of the 

zeroth order inner solution on the initial conditions is expressed in Section 3 through a set 
of 20 equations for 20 unknown parameters that define the inner and outer zeroth order 
solutions. The procedure to obtain the required derivatives involves use of the chain rule, 



where first the partials of /?j(yO with respect to the 20 parameters is taken, and then 
multiplied by the partials of the parameters with respect to the initial conditions. The details 
are omitted, but the procedure is similar to that developed for a regular expansion in Ref. 
11 . 


4.3 First Order Guided Solution 

Gosed loop guided solutions are obtained by using the optimal control expressions 
given in Eq s (3.26) and (3.27). Feedback implementation entails treating the current state 
(from the simulation) at each control update as a new initial state, and calculating the costate 
values corresponding to the same time instant. The estimate for these costates to first order 
is obtained by repetitively following two steps. First the zeroth order MAE problem is 
solved, providing the zeroth order Euler solution which defines the characteristic curves for 
the first order expansion of the HJB equation. The procedure for this step has been 
detailed in Section 3. Next the quadratures in Eq's. (4.63-4.65) are performed to correct 
the costates to first order. When a specified time increment has been reached, the current 
states are used as initial conditions for the next MAE calculation followed by the first order 
correction step. During the exit phase, the left outer solution is discarded. This situation 
represents the case where to>0 in the analysis of Section 4.1. 

Figure 4.2 compares the exact solution, the zeroth order guided solution and the 
first order calculation for the left side (the heading in the entry phase of the maneuver is 
between 0 and 1 1.5°). The exact solution is computed at heading increments of 0.25° 
where the initial conditions are determined by the values of the zeroth order guided 
simulation at the end of the last control update segment. The first order correction is 
calculated but not included in the guided solution. It is clear that initially the first order 
calculation overcorrects the zeroth order solution. Oose to and in the inner region though, 
the first order calculation does provide an excellent correction to the zeroth order solution. 
These partial results indicate the need for further investigation of this behavior. The cause 
may be a numerical ill conditioning in the evaluation of the quadratures in Eq's. (4.63- 
4.65), and perhaps further development of the theory may be required. 






Section V 


Conclusions and Recommendations 


5.1 Conclusions 

In this report it has been demonstrated that the dynamics associated with skip 
trajectories are singularly perturbed, with the perturbation parameter defined as the ratio of 
the atmospheric scale height to a reference radius close to the planet radius. Earlier studies 
have either employed regular perturbation analysis of the inner dynamics, or have 
incorrectly attempted a matched asymptotic expansion (MAE) analysis. We have clearly 
demonstrated that the transformations employed in the earlier regular expansion studies in 
effect transformed the original problem to the inner dynamics, and therefore only have the 
appearance of being the type that occur in regularly perturbed problems. The resulting 
solutions are therefore not uniformly valid. For skip trajectories, this results in a poor 
approximation of the optimal solution near the end of the trajectory, where there is little 
control authority available. 

With regard to MAE analysis, all of the issues improperly dealt with in earlier 
analyses of this type that have been attempted in the past for skip trajectories, have been 
corrected . The first issue deals with the fact that both the inner and outer expansions are 
crucially involved in satisfying the boundary conditions. The use of the outer expansion 
alone to satisfy initial conditions leads to discrepancies in the matching conditions. A 
second issue is the need for separate left and right outer expansions, and the role that the 
inner expansion plays in joining the discontinuities that occur between the outer expansions 
through the matching conditions. The true optimal solution is, by its nature, continuous; 
therefore, in order for the composite solution to serve as a uniformly valid approximation to 
the exact solution, it also must be continuous. In a skip trajectory, a discontinuity between 
the left and right zeroth order outer solutions is caused by the change that occurs in the 
trajectory parameters during the osculating atmospheric maneuver. We have demonstrated 


76 


that the zeroth order inner solution can be viewed as a process whereby the discontinuity 
between the left and right zeroth order outer solutions is taken into account 

A third issue concerns the proper selection of the reference altitude, used in defining 
the independent variable of integration in the outer solution, which has significant 
implications on the solution process. An arbitrary choice of this altitude leads to a situation 
where the outer solution for the flight path angle cannot be evaluated, thereby preventing a 
zeroth order composite solution from being formed. A systematic approach to obtain a 
relationship to determine the correct reference altitude is to use the condition that the 
composite zeroth order flight path angle solution is zero at the lowest point of the trajectory. 

Application of these ideas to the problem of inclination change with minimum 
energy loss has resulted in a zero order solution in the form of a set of 20 algebraic 
equations. By exploiting the structure inherent in the matching procedure it is possible to 
reduce the problem to 6 equations in 6 unknowns. A further simplification was employed 
which permits satisfaction of two of the boundary conditions by partially separating two of 
the equations from the remaining four equations. Numerical experience shows that the 
zeroth order solution is close to the optimal solution, and that the outer solution plays a 
critical role in accounting for the variations in Loh's term near the exit phase of the 
maneuver. This feature is what differentiates a MAE analysis from a regular perturbation 
analysis of skip trajectories. However, the deficiency that remains in several of the critical 
variables indicates the need for a first order correction. 

Expansion of the solution to first order required further development of the MAE 
expansion procedure. We have developed a general procedure for constructing a matched 
asymptotic expansion of the Hamilton-Jacobi-Bellman (HJB) equation, based on the 
method of characteristics. Of particular interest here is the manner in which matching and 
boundary conditions are enforced when the expansion is earned out to first order. We have 
recognized the need to distinguish between two cases pertaining to the location of the 
singular region with respect to the boundary conditions. The first is where the left 
boundary condition coincides with, or lies to the right of, the singular region, and the 
second is where the singular region lies between the boundary conditions. It is shown that 
the boundary conditions for the characteristic curves of the HJB equations are also defined 
by the MAE method. The characteristic curves for the inner and outer expansions of the 
HJB equation are the inner and outer solution components of the zeroth order MAE 
solution. Another consequence of the analysis is that whenever the outer and inner zeroth 
order MAE solutions are used to define the characteristic curves, the first order matching 



conditions are automatically satisfied. A simple example is used to illustrate the procedure, 
where the obtained solution is uniformly valid to 0(e 2 ). The procedure's potential 
application to aeroassisted plane change was also evaluated. 


5.2 Recommendations 

Based on the results developed in this report, the following recommendations are 
made for further research in this area. 

5.2.1 Completion of the First Order Correction Analysis of the 
Aeroassisted Plane Change Problem 

In Section 3 the zeroth order guided solution was obtained and compared with the 
exact solution. The error that exists in some of the zeroth order guided solution variables 
with respect to the exact solution indicated the need for a first order correction. Section 4 
developed a general procedure to obtain a first order correction to the zeroth order problem 
and its potential application to the aeroassisted plane change problem. The numerical 
results indicate the need for further evaluation of the approach, and perhaps further 
development of the theory. Simple examples of increasing complexity proved useful in 
understanding and developing the theory. It is recommended that the first order correction 
for the aeroassisted plane change problem be further analyzed to determine if it is 
numerically ill conditioned, or if an alternative approach or further development of the 
theory is needed. One alternative approach for first order analysis is to expand the Euler 
system equations. 28 


5.2.2 Aerodynamic Healing Requirements 

Aerodynamic heating is an important aspect of the aeroassisted maneuver that was 
not considered in this report To make the guidance law useful in realistic applications, it is 
necessary to include aerodynamic heating requirements in the problem formulation in a 
form such that the zeroth order problem remains tractable. Minimization of the time 
integral of the flight path angle squared may have desirable features, as the numerical work 
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in Ref. 2 indicates. The numerical solution to this problem results in nearly grazing 
trajectories that are considered an useful engineering compromise between energy 
requirements and aerodynamic heating requirements. Hence, a choice of performance 
index that includes a flight path angle squared term is recommended. 

A second approach for formulations that render the zeroth order problem intractable 
(such as enforcing a hard constraint on the heating rate or the heating load) would be to 
investigate combining an analytically tractable portion of the solution with a numerical 
method, such as collocation, to develop a numerically efficient zeroth order solution. The 
approach was developed and illustrated for a launch vehicle application in Ref. 28. 


5.2.3 Atmospheric and Parametric Uncertainties 

The aerodynamic force used to modify the vehicle's trajectory during the 
aeroassisted maneuver is uncertain due to uncertainties associated with estimates of the 
vehicle state vector, atmospheric density and the vehicle's aerodynamic coefficients. The 
effect of these uncertainties along with uncertainties in the entry conditions may result in 
significant trajectory deviations from the nominal trajectory and in large errors in the final 
conditions. 29 An important determinate of the guidance system performance during 
maneuvers, such as aerocapture or landing on the surface of Mars, is the accuracy of the 
information on which the guidance law is based. 3 ^ A further study to improve onboard 
estimation and parameter identification for aeroassisted applications is recommended. 

A second viewpoint is to design the guidance law so that performance is maintained 
in the presence of uncertainty (robust performance). Such a design can be achieved by 
treating the uncertainty as an opponent in a differential game formulation 31 . In the case of 
linear quadratic games, such a formulation is intimately connected to the design of a 
controller that minimizes the infinity norm of the transfer function from disturbances to 
performance outputs. It may also be beneficial to investigate the use of thrust within the 
atmosphere, such as an Aerocruise maneuver to achieve plane change. It has been shown 
in Ref. 32 that Aerocruise is less sensitive to atmospheric uncertainties than aeroglide. 



5.2.4 Expansion of tbe Euler Equations 

An alternative approach to obtaining higher order corrections is to consider an 
expansion of the Euler system of equations. The higher order equations for the states and 
costates are coupled, linear and inhomogeneous, and the contributions of higher order outer 
terms to the composite solution are not zero (as they are in the HJB expansion). The 
solution of these equations requires calculation of a state transition matrix, which can be 
derived analytically using the analytical zeroth order solution, and performing a quadrature 
on both state and costate equations. It has been shown that in the case of a regular 
perturbation, the first order correction obtained by expansion of the HJB equation is equal 
to the correction obtained by expanding the Euler system of equations to first order. 33 
However, this method has an important potential advantage over the HJB expansion in that 
it may be possible to fix the zeroth order solution, and precompute and store the 
quadratures along the zeroth order solution as a function of a monotonic variable (such as 
total energy). Then, at each control update, the state perturbations from the zeroth order 
solution are accounted for in the first order correction by treating them as initial conditions 
for the first order solution. 2 ** This is not possible in an HJB expansion, since the first 
order correction is valid only for the current values of the state and the independent 
variable, and the quadrature must be repeated at each control update. 


5.2.5 Other Problem Formulations 

Extensions of the existing problem formulation would include accounting for the 
effects of the planet’s rotation and cross range angle, which are ignored in the present 
analysis. It is easy to show that these effects are not present in the zeroth order dynamics, 
which means that they appear only in the first order corrections that are computed by 
quadrature. Therefore, conceptually it should be rather straightforward to include these 
effects in the analysis. Also, there are a range of other problems that are of practical 
interest, such as the aerocapture problem and reentry problems. In the aerocapture problem 
the total energy of a vehicle on a hyperbolic trajectory is to be reduced so that it is captured 
by the gravitational field of a planet. Aerodynamic heating limits are an important 
consideration in this problem formulation as well. 
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5.2.6 Parametric Studies 

In the context of parametric studies it is necessary to determine the range of problem 
parameters for which approximate solutions are obtainable, and to conduct a comparison 
with purely exoatmosphenc maneuvers in terms of energy consumption. Typical 
parameters to be considered are the initial and final orbital parameters, and the maximum lift 
to drag ratio. Therefore, it is of interest to conduct evaluations of a more extensive range of 
parameter values than in the limited study conducted here. 



Appendix A 


Relationship for the Zeroth Order Control c 


Expressed in the transformed variables in the inner region, the Hamiltonian function 
has the form given in Eq. (3.55) 

H; = P' 0r 8/c - Pl w f/(j+ p ' 0t ( i + S 2 + (T 2 )/ <r (Al) 

The control functions are obtained from the optimality conditions = 0 and 
H oa = 0- F° r the above expression of HJ, these conditions, after some manipulation, 
are: 


h;, = </t + 2P' 0 ,s)/o = o 


(A2) 


HL = PL - (-Pi 




- Pin + ft )/® 3 = o 


(A3) 


and the control functions are obtained as: 

8 = -P' 0r /2P' 0v (A4) 

02 = (~ P LyI ~ P'o r / 4 P' 0v )/ Pl v + 1 (A5 ) 

Substituting these functions back into the Hamiltonian in Eq. (Al), after some manipulation 
it takes the form: 
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H; = (2PU a)[(-Pl Y l ~ Pj? r l 4Pl)/Po, + 1] = 2 aP' Qv (A6) 

where the expression in the brackets is recognized to be o 2 . Since HJ, and P' 0v are 
constant in this problem, so is a. Thus it should be possible to express a in term of 

constant parameters only, which is not immediately recognized in Eq. (A5) since it contains 
the functions P' 0r and y ' 0 . To show this, use of Eq. (3.64) in Eq. (A5) yields 

+ cyiPlf - pIyjpI + l (A7) 

From Eq. (3.67) P' 0w and c can be written as: 

PL = (A8) 

c = -2 aPlk 2 (A9) 

Using these equations and Eq. (3.61) in Eq. (A7), a becomes: 


a = 1/(1 + k \ + 2*,* 3 ) U2 


(A10) 


Appendix B 


Initial Guess and Evaluation 
of the Inner Integration Constants 


B.l Initial Guess Procedure 


To calculate the 6 inner integration constants k 3, k4, k5, c, P' 0w , and P' 0v in Eq's. 

(3.61-3.66), 6 boundary values must be supplied. These are the left common part values 
of y, v and \y and the right common part values of y, P v and y which will be denoted by 

the subscripts L and R" respectively. At the very first step of the procedure, the actual 
left boundary values of y, v and \jr and the right boundary values of y and y are used as an 
initial guess for y L , v L , y/ L , y R and yr R respectively. The value of P v , is estimated using 
the boundary condition in Eq. (3.1 16) that /*(*,)= 1.0. Since h f is far from the region 
of singularity (h= 0), this is closely approximated by the outer solution, so R P® U (hj- ) = 1.0. 
Since h« 1, the value of (h) does not change substantially in Eq. (3.31), so 
w o (hf) an ^ from Eq. (3.38) R Po«(0)= R / > ) J| l (/iy) = 1.0. These properties can be 

verified by examining the converged solution for P° u (y/) in Fig. 3.7. Next the 
transformation in Eq. (3.83) is used to relate R P 0 ° U ( 0) to R P 0 ° V (0), where R « 0 o (0)= R u 0 °(/ I/ ) 
and from Eq. (3.44) 


V(A/> 


= < r Rv o° (*/)/£* 


I R 


s r s 


(Bl) 


Now v 0 (hf)~ vl(\f/jr)- v„(^y) and Vq(^) is calculated by approximating the 
dynamics in Eq. (3.49) using dv\ / dy = 2.0. All of the above properties may be verified 
from Fig. 3.5. Finally, the matching condition in Eq. (3.86)P 0 ° V (0) = pj v ( ~) = P VR i s used 
to provide an estimate for P VR . 
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B.2 Inner Integration Constants 

Using Eq's. (3.61) and (3.62) at the left and right common part boundary 
conditions, kj, k 2 and k 3 become: 

= -12oAw/Ayr 3 - 6 (y R + y^/Ayr 2 (B2) 

*2 = Ay/ Ay/ + *,(y r R + ^)/2 (B3) 

* 3 = ~KWrWl / 2 + (YlVr ~ YrVi)/AY (B4) 

where A( ) = (•)* - ( ) t and a is found as a function of the boundary conditions by 
using Eq’s. (B2-B4) in Eq. (A 10). The expression for a becomes 

<7 = -^Ayr 2 /[Ayr 2 +3(y R + y L ) 2 + Ay 2 ] (B5) 

From Eq. (3.43) at 77 — > 00 , the common part values of w are zero, thus 
w L = w R =0 when the initial conditions are in the left side (entry phase). In this case Aw=0 
in Eq. (B.2). When the initial conditions are given in the right side (exit Phase), the left 
matching condition is dropped and vv L is equal to the current w. Use of Eq's. (B2) and 
(B3) and the estimate of P[ in Eq. (3.67) provides the values of the constants P' 0w and c. 

Next, k 4 and k 5 can be found using Eq's. (3.62) and (3.63) at the left boundary 
conditions: 

k 4 = ow L — k^y/ L / 6 + k^yfiH + k^y/ L (B6) 

*5 = v t - (0 + l/< j)y l - o[(yr t *, ~ * 2 ) 3 ]/ 3 ^, (B7) 

At this stage all the inner states and costates can be evaluated at any y between \j/l 
and \|Tr as a function of the boundary conditions using Eq's. (3.61-3.66). In particular, the 
initial and final inner solutions are needed to construct the composite solution, which is 
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used to enforce the boundary conditions. Also, at the lowest point in the trajectory (h=0, 
r=r s ) the composite solution is equal to the inner solution and the flight path angle yj(0,e) 

is zero. Thus, the value of \|/ s that corresponds to this point can be found from Eq. (3.61 ) 
by equating it with zero, and then using the result in Eq. (3.62) to calculate the 
corresponding value of w s . This point is illustrated in Fig. 3.2. Finally, p s and hence, the 
reference radius r s , can be obtained by using w s and the relationships defined in Eq. 
(3.43). 


V 


V 
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Appendix C 

Relationship Between the Zeroth 
and the First Order Matching Conditions 


The quantity R ,'(r) needed for the quadratures in Eq. (4.59) is given by Eq. (4.15) 
R! = -PL(T)fti*( r).0) (Cl) 

where x(t) satisfies Eq. (4.17). Therefore 

R'r = -/t(~)4(4(~),o) (C 2 ) 

where x' 0 (z) and P' 0 x (t) are the state and costate solutions along the inner characteristic 
curve. The right side matching condition in Eq. (4.25) imposes a constraint that must be 
satisfied by the inner and outer characteristic curves. That is 

/l(~)/oW~)>0) = *C(0)/o(%°(OXO) (C3) 

where the superscript R is used in Eq. (C3) in recognition of the fact that the outer 
characteristic curves are discontinuous at t=0. Noting that f' 0 and / 0 ° are the same 

functions and use of Eq. (C3) together with the matching conditions in Eq. (4.23) imply 
that: 


RU~) = ^ 0 ,( 0 ) . 


4(oo) = X(0) 


(C4) 
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V-' 


A similar consideration of the left side matching conditions in Eq's. (4.54) and (4.56) leads 
to the conditions: 


= g K<P) 


x { 0 (-~) = l *o(0) 


(C5) 


It can be seen that the inner and outer characteristic curves must also satisfy the MAE 
matching conditions of Eq. (2.5). Thus it follows that the Euler solutions from a zeroth 
order MAE analysis, as outlined in Section 2, serve as the characteristic curves along which 
the quadratures in Eq. (4.59) are performed. 


W 
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